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ABSTRACT 

A   simple  physical    model    of   a    copper    crystal    surface   was    de- 
veloped.     Atoms    were   considered  as    quasi-hard   spheres    which 
occupied  perfect    lattice  positions.       A   computer    simulation,    based 
on    energy   consideration   only,    using    the   Monte   Carlo  method   was 
developed,    tested  and  used   to   study    equilibrium   surface   micro- 
states.      As    a    result    of    this    study,    four    conclusions    were    drawn: 

1.  This   model    holds   promise    for    further    investigation    of    real 
crystal    surface  phenomenon. 

2.  Minimum   energy    considerations    cause   atoms    to  align    them- 
selves  preferentially    in   a    (lio)    direction    on    the   surface    of  a 
face-centered  cubic    crystal. 

3-       Stepped   surface   configurations    are    fairly    stable,    but 
isolated   "stub'^   atoms   and   vacancies    tend   to   coalesce   with    other 
""stubs"   and   vacancies,    respectively. 

4.      Random  motions    of   the    individual   atoms    cause   aggregates    of 
atoms    to  break  apart   and  recombine. 
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I.     I^r^RODUCTION 

A.   HISTORY 

This  simulation  was  based  on  two  areas  of  contemporary  research; 
crystal  growth  and  radiation  damage.   Ideas  were  borrowed  from  each 
of  these  fields  of  study  in  order  to  devise  a  new,  simple  physical 
model  of  a  real  crystal  surface,  and  to  begin  the  study  of  crystal 
surface  equilibrium  microstates.   Brief  discussions  of  prior  re- 
search is  contained  in  the  following  paragraphs. 

Theoretical  models  of  crystal  growth  have  been  studied  for  over 
20  years,  and  numerous  papers  have  been  written  attempting  to  de- 
velop models  for  growth  from  a  solution,  a  melt,  or  a  vapor 
[4,6,9,10].        An  essential  part  of  these  investigations  was  the 
assumption  of  a  perfect  simple  cubic  crystal  lattice  and  of  first 
nearest  neighbor  forces  only.   Some  results  have  been  compared  with 
experimental  work  done  on  the  growth  and  dissolution  of  simple 
cubic  ionic  crystals  in  supersaturated  and  unsaturated  solutions, 
and  have  shown  both  agreement  and  disagreement  with  the  theory  [6] . 
Motion  of  atoms  absorbed  on  the  surface  of  the  crystal  were  taken 
into  account  by  incorporation  of  theoretical  surface  diffusion  mech- 
anisms based  on  transport  theory.   These  papers  were  generally 
written  in  terms  of  classical  thermodynamic  quantities.   Crystal 
growth,  diffusion  and  dissolution  were  all  discussed  in  terms  of 
Gibbs  free  energy  and  chemical  potential  for  generalization  to  any 
simple  cubic  crystal.   The  crystals  visualized  were  usually  ionic. 


Radiation  damage  studies  provided  useful  physical  information 
on  real  crystals.   Gibson,  Goland,  Milgram  and  Vineyard  (hereafter, 
GGMV) ,  and  Girifalco  and  Weizer  (hereafter,  GW)  have  investigated 
interatomic  potentials  which  were  found  useful.   GGMV  used  a  Born- 
Mayer  exponential  approximation  to  the  interatomic  potential,  and 
applied  their  results  to  radiation  damage  dynamics  for  short  in- 
teratomic spacings  CSl •   The  GGMV  potential  used  here  has  become 
known  as  the  "Gibson  II  potential".   GW  used  a  Morse  potential  and 
fitted  data  for  15  cubic  metals  to  arrive  at  a  appropriate  Morse 
potential  parameters  for  metals  [?] •   Anderman  followed  the  GW 
procedure  and  obtained  slightly  different  parameters  for  use  in 
the  Morse  potential,  by  introducing  the  idea  of  truncating  the 
potential  function.   Anderman  used  his  values  in  some  static  radia- 
tion damage  simulations  [l].   The  Anderman  potential  values  have 
been  used  here,  and  the  Anderman  form  of  a  GW  Morse  potential  will 
be  referred  to  as  a  "truncated  Morse  potential". 

Since  solving  for  the  motion  of  the  atom  in  a  crystal  involves 
solving  an  N-body  problem,  an  essential  assumption  for  the  radia- 
tion damage  studies  was  that  the  motion  of  an  atom  could  be  des- 
cribed by  an  independent  set  of  two-body  potentials  such  as  the 
Gibson  II  and  truncated  Morse  potentials.   Such  a  technique  yielded 
a  soluble  approxmation  to  the  actual  physical  problem.   For  ex- 
ample, an  atom  and  four  of  its  first  nearest  neighbors  gave  four 
independent  two-body  problems,  instead  of  one  five-body  problem. 

B.   THE  NEED  FOR  NUMERICAL  SIMULATION 

The  detailed  processes  involved  in  both  of  these  areas  of  study 
were  essentially  unknown.   Both  lines  of  research  were  attempts  to 
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take  macroscopic  information  and  develop  microscopic  models  which 
predict  macroscopic  phenomenon.   In  some  cases  the  theory  involved 
did  not  permit  analytical  solution  of  the  problem  [6] 
and  in  all  cases  random  processes  at  the  atomic  level  were  assumed. 
Solving  the  same  problem  a  number  of  times  using  a  Monte  Carlo  model 
and  averaging  the  properties  of  several  final  configurations  was  a 
technique  common  to  both  of  these  areas  of  interest.   Since  the 
calculations  were  repetitive  and  could  be  programmed  fairly  simply, 
computer  simulation  was  a  logical  tool  to  use.   In  addition,  since 
microscopic  models  were  assumed  and  system  microstates  were  gene- 
rated, a  computer  simulation  yielded  microscopic  information  - 
which  may  be  appropriately  averaged  and  interpreted  to  yield  theo- 
retical macroscopic  results  for  comparison  with  empirical  macro- 
scopic results.   A  "computer  experiment"  therefore  yielded 
microscopic  as  well  as  macroscopic  information. 

C.   THE  MONTE  CARLO  METHOD 

The  Monte  Carlo  method  requires  that  in  calculation,  whenever 
one  makes  a  decision  based  on  probability,  the  decision  is  made  by 
generating  a  random  number  and  comparing  it  with  the  probability 
involved.   The  calculation  is  continued  until  it  arrives  at  another 
decision  based  on  probability.   The  method  then  requires  the  gene- 
ration of  a  new  random  number  for  comparison  and  decision.   The 
entire  process  continues  for  as  long  as  desired,  and  a  set  of  final 
conditions  is  obtained  from  a  set  of  given  initial  conditions.   The 
entire  process  may  then  be  repeated,  starting  with  the  original  set 
of  initial  conditions  and  new  random  numbers,  to  arrive  at  a  new 
set  of  final  conditions.   In  the  end  there  is  an  ensemble  of  final 
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states  generated  from  a  single  initial  state,  and  it  is  possible 
to  compute  the  average  properties  of  the  final  state  L12].   The 
average  final  state  properties  may  be  compared  with  experimental 
data  to  yield  information  on  the  quality  of  the  theoretical  micro- 
scopic model.   Essentially,  the  method  attempts  to  generate  macro- 
scopic information  from  a  microscopic  system  model. 

D.   THE  MODEL  ATTEMPTED  HERE 

The  present  effort  is  believed  to  be  the  first  attempt  to  study 
the  surface  dynamics  of  an  actual  metal  crystal  affected  only  by 
the  temperature  of  the  system.   No  attempt  was  made  to  include 
growth  or  dissolution  (adding  or  subtracting  atoms)  of  a  crystal. 
Rather,  a  simple  physical  approach  was  tried  in  an  effort  to  study 
only  the  surface  of  an  actual  metallic  crystal.   What  was  visualized 
was  a  section  of  a  perfect  copper  crystal  surface  that  was  initially 
"piled"  in  different  ways.   Observations  were  made  of  the  subsequent 
surface  configurations  (micr estates )  resulting  from  Monte  Carlo 
calculations.   Interatomic  potentials  were  borrowed  from  the  radi- 
ation damage  studies.   Transition  probabilities  were  borrowed  from 
the  thermodynamic  crystal  growth  studies.   An  appropriate  computer 
program  was  written  to  simulate  the  copper  surface  dynamics,  and 
several  computer  runs  were  made  to  obtain  a  rough  idea  of  what 
appeared  to  be  equilibrium  microstates.   The  programming  language 
used  throughout  this  study  was  FORTRAN  IV,  and  the  computer  used 
was  the  IBM  36O/67. 
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II.   THE  SIMULATION  MODEL 

A.   ASSUMPTIONS 

Since  this  was  a  first  attempt  at  solving  a  new  problem, 
several  simplifying  assumptions  were  made.   Perhaps  the  greatest 
assumption  made  in  this  calculation  was  that  metal  crystal  binding 
energies  could  be  computed  from  a  single  composite  potential  energy 
function,  with  the  interactions  between  atoms  considered  as  in- 
dependent two-body  problems.   The  essential  feature  of  this  as- 
sumption was  to  treat  the  atoms  of  a  crystal  lattice  as  quasi-hard 
spheres  which  interacted  by  some  "unknown"  mechanisms  to  yield 
observed  crystal  behavior. 

The  only  metal  visualized  was  copper ,  which  forms  a  face- 
centered  cubic  (fee)  lattice.   For  simplicity,  the  (001)  crystal 
orientation  was  assumed.   Lattice  sites  were  visualized  as  ordered 
real  triples  in  the  positive  octant  of  a  rectangular  coordinate 
system.   The  coordinate  planes  and  axes  bordering  the  positive 
octant  were  considered  as  belonging  to  the  positive  octant.   See 
Figure  1.   Although  a  rectangular  parallelepiped  was  always  gene- 
rated to  provide  possible  lattice  locations  (hereafter  called  the 
"active  lattice  volume"),  appropriate  changes  in  the  lattice  gene- 
rator insured  that  the  initial  configuration  of  the  surface  could 
take  any  desired  shape  with  the  sole  precondition  that  the  only 
sites  available  were  perfect  lattice  sites.   The  active  lattice 
volume  was  considered  as  merely  a  framework  of  sites  that  were 
available  for  atoms- whether  the  sites  were  actually  occupied  or  not 
depended  on  the  details  of  the  microstate. 
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In  addition,  it  was  assumed  that  only  the  top  several  atomic 
layers  of  the  surface  provide  interesting  results  -  in  effect,  this 
meant  that  the  active  lattice  sites  generated  were  placed  on  top  of 
a  perfect  lattice  substructure  which  was  not  allowed  to  interact 
with  the  surface  sites  in  any  way  other  than  to  be  used  in  calcu- 
lating surface  site  potential  energy.   This  further  implied  that 
above  the  active  lattice  sites  all  locations  were  always  kept 
vacant . 

Several  additional  assumptions  were  made  to  simplify  the  model. 
These  were: 

1.  Atomic  vibrations  and  zero-point  energy  were  neglected. 

2.  Lattice  relaxation  due  to  vacancies  was  neglected. 

3.  No  interstitial  atoms  were  allowed. 

4.  No  impurity  atoms  were  allowed. 

5.  Only  even  numbers  of  planes  in  the  principal  coordinate 
directions  were  allowed. 

Further  discussion  of  these  assumptions  is  contained  in  Appendix  A. 

B.   THE  COMPUTER  PROGRAM 

A  computer  program  that  provided  an  adequate  picture  of  micro- 
scopic surface  activity  had  to  be  written  and  tested.   Overall 
program  flow  is  shown  in  Figure  2.   The  main  features  of  the  pro- 
gram are  discussed  in  the  remainder  of  this  chapter.   Further 
details  are  provided  in  Appendix  A. 

1.   The  Lattice  Generator 

The  first  task  involved  generating  an  appropriate  lattice- 
work of  crystal  sites  in  the  computer.   For  the  (001)  orientation 
of  a  face-centered  cubic  crystal  (hereafter,  "fee,  (001)  crystal"), 
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available  sites  are  described  by  coordinates  such  that  the  sum  of 
coordinates  is  either  an  even  or  an  odd  integer  for  all  positions 
in  the  crystal.   The  even  choice  was  used  in  this  simulation. 

After  a  valid  lattice  site  was  located,  it  was  given  a 
number  as  a  label.   Then  an  occupation  index  was  defined  for  the 
site  to  retain  a  knowledge  of  the  surface  configuration.   Several 
different  initial  configurations  were  used,  and  these  are  sketched 
as  the  "a"  parts  of  Figures  3  through  1?.   The  specific  program 
statements  to  generate  such  surfaces  are  contained  in  Appendix  B. 
2 .   Periodic  Boundary  Conditions 

Computer  storage  and  time  limitations  made  impractical  the 
use  of  an  arbitrarily  large  active  lattice  volume.   A  practical 
lattice  size  involved  at  most  a  few  thousand  active  lattice  sites. 
This  required  that  the  active  volume  size  be  limited  to  a 
Cartesian  space  of  roughly  20  x  20  x  10,  which  corresponds  to 
2,000  active  lattice  sites  and  a  copper  crystal  measuring  about 
50A  x  50A  X  25A.   Such  a  small  crystal  with  six  surfaces  compli- 
cated the  study  of  surface  effects. 

One  surface  (the  bottom  of  the  active  volume)  has  already 
been  eliminated  by  assuming  that  the  active  lattice  volume  was 
superimposed  on  a  perfect  crystal  substrate  which  was  never  allowed 
to  vary.   Four  other  surfaces  (the  sides  of  the  active  volume)  were 
eliminated  by  assuming  "periodic  boundary  conditions"  which  resulted 
in  an  essentially  infinite  single  surface  (the  "top"  of  the  active 
volume)  that  varied  periodically  in  configuration.   This  infinite 
surface  may  be  visualized  as  part  of  a  crystalline  surface  far 
from  a  grain  boundary. 
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These  periodic  boundary  conditions  were  achieved  by 
visualizing  that  each  active  volume  plane  was  laterally  surrounded 
by  exactly  similar  active  volume  planes  that  had  been  translated  in 
the  X  and/or  y  dir ect ion( s ) .    In  other  words,  an  atom  on  the  left 
edge  of  an  active  volume  plane  "sees"  the  right  edge  of  that  same 
active  volume  plane;  an  atom  on  the  top  edge  "sees"  the  bottom  edge 
of  the  same  plane,  and  similarly  for  the  other  edges. 

Periodic  boundary  conditions  were  achieved  in  practice  by 
defining  nearest  neighbor  (NN)  arrays  for  each  site.   These  arrays 
were  calculated  by  visualizing  the  spatial  relationships  between 
an  atom  in  the  fee,  (001)  crystal  and  its  nearby  surrounding  sites - 
see  Figures  l8  and  19.   Atoms  within  one  atomic  layer  of  the  sur- 
face of  the  active  lattice  volume  required  special  attention  to 
insure  the  proper  periodic  structure  of  each  plane. 
3.   Potential  Energy 

The  model  chosen  assumed  that  the  crystal  may  be  approxi- 
imated  by  a  latticework  of  quasi-hard  spheres,  so  a  scheme  for 
holding  individual  atoms  together  was  needed.   The  problem  was  an 
N-body  problem,  and  therefore  no  general  exact  solution  existed. 
However,  approximating  the  potentials  by  an  independent  set  of  two- 
body  potentials  yielded  a  solvable  approximation  to  the  physical 
system,  and  this  was  a  standard  technique  used  in  radiation  damage 
studies . 

The  next  task  involved  a  judicious  choice  of  potential 
function.   A  composite  potential  function  was  chosen.   For  low 
values  of  interatomic  spacing  the  Gibson  II  potential  Csl  was  used, 
and  for  larger  values  the  truncated  Morse  potential  Cl]  was  used. 
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since  these  two  potentials  do  not  match  smoothly  together  (see 
Figure  20),  a  cubic  function  was  generated  which  matched  the  value 
and  the  slope  of  these  two  potentials  at  appropriate  points  so 
that  a  smooth  composite  potential  was  achieved  (see  Figure  21). 
The  cubic  potential  coefficients  were  generated  by  solving  a  system 
of  four  simultaneous  linear  equations  using  Subroutine  CROSYM 
developed  at  Lockheed.   The  form  of  the  potential  equations,  their 
appropriate  coefficients,  and  the  cutoff  values  used  in  this  simu- 
lation are  listed  in  Table  I. 

Next,  a  potential  energy  was  associated  with  each  site  in 
the  lattice  based  on  the  number  of  its  first  and  second  nearest 
neighbor  (NNl  and  NN2 ,  respectively)  sites  that  were  occupied.   In 
other  words,  every  active  volume  lattice  site  was  given  a  potential 
energy.   If  the  site  were  occupied,  this  energy  was  just  the  total 
mutual  interaction  energy  between  an  atom  and  its  NNl's  and  NN2's. 
However,  if  the  site  were  vacant,  this  energy  was  the  interaction 
energy  between  the  NNl's  and  NN2 '  s  of  the  site  that  would  arise  i_f 
the  site  were  occupied.   In  practice,  potentials  were  given  to  each 
site  (occupied  or  vacant)  in  the  same  manner,  but  the  use  to  which 
these  potentials  were  put  depended  vitally  on  whether  a  particular 
site  was  occupied  or  not. 

On  the  average  about  half  of  the  mutual  interaction  energy 
between  sites  will  be  associated  with  each  site  individually,  so  a 
parameter  PFIV  (mnemonic  for  "point  five",  and  usually  equal  to 
0.5)  was  used  as  an  energy  distribution  factor.   By  varying  this 
energy  distribution  factor,  physical  effects  resulted  which  were 
equivalent  to  varying  the  lattice  temperature.   Therefore,  by 
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changing  the  energy  distribution  factor,  different  lattice  tem- 
peratures were  achieved.   The  relationship  between  the  equivalent 
lattice  temperature  and  PFIV  is  indicated  in  Figure  22  and  dis- 
cussed in  greater  detail  in  Appendix  A. 
4.   Transition  Probability 

If  a  given  lattice  site  were  occupied  and  had  one  or  more 
vacant  sites  as  NNl's  it  had  (statistically)  a  finite  probability 
of  moving  from  its  present  site  to  the  nearby  vacant  site.   In 
practice,  when  an  occupied  lattice  site  was  reached,  the  computer 
searched  for  lower  values  of  site  potential  among  the  unoccupied 
NNl's  of  the  occupied  site  of  interest.   Such  lower  energy  vacant 
sites  implied  that  the  atom  under  consideration  was  sitting  on  the 
edge  of  a  potential  well,  and  a  transition  was  carried  out  to  the 
deepest  NNl  well  -  that  is,  the  atom  "jumped"  into  the  deepest 
NNl  well.   If  a  "deepest  well"  was  not  unique,  then  the  actual 
final  site  was  chosen  stochastically.   "Wells"  of  zero  depth  were 
considered  as  having  "lower  energy"  for  the  purpose  of  this 
calculation . 

For  atoms  with  no  vacant  NNl  sites  of  lower  energy  but 
one  or  more  with  higher  energy,  the  atom  would  have  to  climb  a 
potential  hill  to  reach  it  (them).   This  was  considered  possible 
but  was  given  a  probability  defined  by  a  Boltzmann  factor: 

p  =  exp  (-E/kT)  (1) 

where 

P  =  the  probability  of  transition  to  a  higher  energy  site 

£(-^0)  =  energy  difference  between  the  occupied  site  of 

interest  and  one  of  its  NNl  sites  with  higher  energy 
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k  =  Boltzmann's  constant 
T  =  absolute  temperature 

The  computer  searched  for  the  lowest  positive  energy  difference 
and  made  the  decision  to  jump  or  not  stochastically.   If  several 
NNl  sites  had  positive  energy  differences,  they  were  considered 
in  turn  from  the  lowest  to  the  highest  until  a  transition  was  made. 
If  no  transition  were  made  by  the  time  the  positive  energy  dif- 
ferences were  exhausted,  then  the  atom  stayed  where  it  was  for 
that  particular  microstate. 

After  any  transition  the  potential  energy  of  the  NNl's  and 
NN2's  of  both  the  old  and  new  sites  were  adjusted  to  account  for 
the  changes  made.   After  every  third  microstate,  all  potentials 
were  recalculated  to  avoid  loss  of  accuracy.   A  pseudo-random 
number  generator  was  used  to  produce  the  random  probabilities 
needed  to  make  the  decisions  involved. 
5.   Computer  Output 

The  problem  of  how  to  present  microstate  information  in 
easily  interpreted  form  was  solved  by  having  the  computer  "draw" 
pictures  of  microcrystal  planes,  called  "microstate  pictures". 
The  computer  printed  arrays  of  occupation  indices  to  indicate  a 
vacancy  or  an  occupied  site,  and  each  digit  was  in  its  correct  x 
and  y  coordinate  position  for  the  site  it  represented.   See 
Figure  23. 
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III.   RESULTS 

A.   INITIAL  TESTING 

Test  runs  made  prior  to  taking  any  "good  data"  indicated  that 
two  problems  would  be  encountered.   The  first  problem  discovered 
was  the  lack  of  net  motion  out  of  a  perfect  surface  configuration. 
The  second  problem  was  an  indication  that  the  model  might  tend  to 
"transport"  atoms  preferentially  in  any  given  plane  from  higher  to 
lower  values  of  the  arbitrarily  assigned  y-coor dina te. 

The  problem  of  no  net  motion  from  a  perfect  surface  was  solved 
by  two  different  techniques.   One  was  to  alter  the  energy  distri- 
bution factor  PFIV  to  define  an  equivalent  lattice  temperature  to 
literally  force  some  action.   The  other  method  involved  piling  the 
surface  atoms  in  both  random  and  preconceived  configurations  for 
the  initial  microstate  and  observing  subsequent  microstates  using 
the  correct  average  value  of  PFIV  (=  0.5).   The  results  of  these 
tests  are  discussed  in  Sections  B  and  C,  below. 

The  problem  of  model  dependent  transport  was  more  fundamental, 
since  any  results  obtained  would  be  biased  in  some  unknown  fashion, 
Again,  two  techniques  were  used  to  explore  the  problem.   The  first 
technique  was  to  test  the  motion  of  a  single  atom  initially  placed 
on  a  perfect  surface.   The  motion  of  this  single  "stub"  atom  was 
observed  for  400  microstates.   The  motion  within  each  microstate 
was  recorded  as  the  net  number  of  lattice  units  moved  "up",  "down", 
"left"  or  "right"  from  its  previous  position.   These  directions 
were  defined  for  an  observer  looking  at  a  single  lattice  plane; 
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"right"  meant  "in  the  +  x-dir ect ion " ,  and  "up"  meant  "in  the 
+  y-direction"  for  the  plane  under  consideration.   The  average 
results  are  reported  in  Table  II  to  the  nearest  one  percent.   A 
quick  glance  at  the  figures  appeared  to  lead  to  the  conclusion  that 
the  model  was  biased  to  motion  right  and  down.   However,  consider- 
ation of  the  average  number  of  jumps  weighted  by  the  distance 
jumped  indicated  that: 

1.  When  the  stub  moved  left,  its  average  motion  was  0.47 
lattice  units . 

2.  When  the  stub  moved  right,  its  average  motion  was  0.6l 
lattice  units . 

3.  When  the  stub  moved  up,  its  average  motion  was  0.46 
lattice  units . 

4.  When  the  stub  moved  down,  its  average  motion  was  0.44 
lattice  units . 

This  meant  that  the  average  net  motion  of  a  single  stub  atom  was 
right  and  almost  entirely  lateral. 

Since  this  test  was  performed  with  only  one  stub  atom,  a  second 
technique  was  used  to  test  for  model  biasing.   The  technique  was 
to  study  the  apparent  motion,  if  any,  of  large  numbers  of  atoms 
placed  on  a  perfect  surface.   This  required,  in  effect,  "taking 
data"  and  the  results  were  available  only  after  several  long  com- 
puter runs  were  made.   Based  on  the  runs  discussed  below,  it  was 
concluded  that  there  existed  no  significant  motion  of  large  numbers 
of  atoms  attributable  to  model  biasing.   Some  motion  "down"  was 
found,  but  this  amounted  to  about  one  lattice  unit  every  10  to  20 
microstates  for  100  or  more  atoms  and  was  not  considered 
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significant.   This  observation  was  purely  qualitative,  and  such 
motion  was  not  found  in  every  case.   In  addition,  no  gross  lateral 
motion  was  observed,  and  any  that  did  exist  was  much  slower  than 
the  "down"  motion. 

B.   PERFECT  SURFACES 

Tests  made  using  perfect  surfaces  required  that  the  energy 

distribution  factor  PFIV  be  altered.   The  arbitrary  initial  tem- 

o 
oerature  chosen  for  all  calculations  was  1,000   K.   With  the 

correct  average  value  of  PFIV  (=  0.5),  no  net  change  was  observed 
for  any  microstate  in  500  attempted  microstates  when  the  initial 
microstate  was  a  perfect  surface. 

In  an  effort  to  force  some  different  microstates  to  occur, 
PFIV  was  lowered  by  a  factor  of  10,  yielding  an  equivalent  tem- 
perature of  10,000   K  -  well  above  the  copper  boiling  point  of 
2,583   K.   This  resulted  in  a  decidedly  unphysical  sequence  of 
surface  microstates.   Since  the  temperature  was  so  high,  the  sur- 
face literally  evaporated  and  burroughed  a  deep  hole  into  itself 
within  10  microstates.   The  model  restraints  insured  that  no  atom 
would  be  lost,  but  this  did  not  mean  that  the  microcrystal  did  not 
try  to  "explode"  -  a  program  flag  included  to  warn  that  an  atom  had 
no  NNl's  was  flashed  at  least  35  times  during    simulation  of  500 
microstates  and  never  appeared  in  any  other  simulation.   Equilibrium 
numbers  of  atoms  in  each  plane  was  reached  after  about  3O  micro- 
states  (as  opposed  to  more  than  100  microstates  normally  required 
for  other  simulations).   Some  single  atom  vacancies  and  stubs  were 
observed,  but  such  surface  defects  were  repaired  quickly.   However, 
the  frequency  of  observation  of  these  defects  was  much  greater  in 
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this    simulation    than    for    all    other    simulations.      Figure   24    shows 
a    typical    microstate   achieved.      Note   that    the   figure    is    drawn   with 
lines    bordering    the   atoms   and   atom   clusters    in   a    (lOO)    direction 
(Figure   24a)    and    in   a    (llO>    direction    (Figure   24b).       For    all    other 
simulations,    lines    parallel    to   a    (lio)    direction   appeared   to   be 
most    natural    to   a    face-centered  cubic    crystal,    but    here   the   jagged 
edges   appearing    in   both   sketches    indicated   that    neither    was 
"natural".       Chapter    IV  contains    a    discussion   of   the    "natural" 
crystal    orientation.         Also   note    that    there   are   several    "overhangs" 
visible    in    each   sketch    -   overhangs    were    observed    in    other    simu- 
lations,   but    not    with   the   great    frequency   with   which    they   appeared 
here. 

The   next   perfect    surface  attempted  was    at    an   equivalent    tem- 
perature   of    1,355     K    -    just    one    degree  Kelvin   below   the    copper 
melting  point.       In   200    raicrostates    no  net    change  was    observed   for 
any   microstate. 

Next,    a  perfect   surface   at   an   equivalent    temperature    of 
2,000      K   was    attempted.      A    typical   microstate   achieved   is    indicated 
in   Figure   25.      The  atoms   making   up    the    top   perfect    surface    layer 
"jumped  up"    and    coalesced    on    the   surface. 

Finally,    a    simulation   at    2,584      K    (one    degree   Kelvin   above    the 
copper    boiling  point)    was   attempted.       Equilibrium   numbers    of   atoms 
in   each   plane  were    reached  by   about    microstate    100.      A   typical 
microstate    is    indicated    in    Figure   26.      Atoms    in   several    surface 
layers    jumped  up   and  actually    formed  steps    with   edges   parallel    to 
the    (no)    direction. 
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C.   PERFECT  SURFACES  PLUS  EXTRA  ATOMS 

Most  simulations  attempted  fell  into  this  category.   This  type 
of  initial  configuration  was  used  for  two  reasons: 

1.  It  was  unphysical  to  expect  that  real  crystal  surfaces 
would  be  absolutely  perfect  for  large  sections  of  surface. 

2.  It  was  unphysical  to  require  temperatures  greater  than  the 
melting  point  in  order  to  achieve  some  idea  of  what  surface  equili- 
brium raicrostates  were. 

All  of  these  simulations  were  run  at  a  lattice  temperature  of 
1,000  K  and  using  the  correct  average  value  of  the  energy  distri- 
bution, PFIV  (=  0.5) . 

Two  simulations  were  run  placing  atoms  randomly  on  top  of  a 
perfect  surface.   Different  random  number  seeds  were  used  for  each 
simulation,  resulting  in  different  numbers  of  extra  atoms  on  top 
of  a  perfect  crystal  surface.   See  Figures  3a  and  4a.   The  randomly 
placed  atoms  coalesced  within  20  micr estates  in  both  simulations. 
Figures  3b  and  4b  show  typical  equilibrium  configurations  for  each 
experiment.   Note  that  the  atoms  arranged  themselves  preferentially 
along  the  (llO)  direction. 

One  simulation  was  run  placing  an  extra  half  plane  of  atoms  on 
top  of  a  perfect  surface,  with  the  halfplane  edge  parallel  to  the 
X-axis.   See  Figure  5a.   At  equilibrium,  the  atoms  were  still  to- 
gether but  again  showed  a  preference  for  lining  up  parallel  to  the 
(no)  direction.   See  Figure  5b. 

Two  simulations  were  run  piling  atoms  into  a  pyramid  config- 
uration with  sides  of  different  slopes.  Sketches  of  the  pyramid 
initial  microstates  may  be  found  as  Figures  6a  and  7a  with  the 
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atoms  aligned  in  the  (lOO)  and  (OlO)  directions.   The  final  con- 
figurations are  shown  in  Figures  6b  and  7b.   Note  that  in  Figure  6b 
all  piled  up  atoms  have  fallen  down  to  coalesce  on  top  of  the 
original  base  plane.   Note  that  in  Figure  7b  not  all  of  the  atoms 
have  fallen  down  -  in  fact,  there  are  aggregates  of  atoms  in  four 
planes  above  the  original  base  plane.   The  large  "hole"  appearing 
in  the  plane  above  the  base  plane  is  the  outer  "groove"  appearing 
in  the  original  microstate  that  has  been  reoriented  and  partially 
filled  in.   Note  that  one  atom  has  actually  jumped  out  of  the  base 
plane  and  has  not  been  replaced  -  this  is  represented  by  the  "hole" 
in  the  left  corner  of  Figure  7b. 

Two  simulations  were  run  with  the  atoms  arranged  to  form  a 
"pyramid  hole"  of  different  step  widths  -  see  Figures  8a  and  9a. 
The  configurations  500  microstates  later  are  shown  in  Figures  8b 
and  9b.   Note  that  in  both  cases  the  atoms  have  fallen  down  into 
lower  planes,  but  that  the  "falling  down"  process  is  a  bit  slower 
when  the  hole  step  is  2  atomic  layers  wide  then  it  is  for  a  hole 
step  of  3  atomic  layers.   The  reason  for  this  is  probably  the  fact 
that  with  2  atomic  layers  per  step  there  are  more  atoms  in  the 
active  lattice  volume,  and  the  active  volume  itself  is  larger. 

Three  simulations  were  run  piling  atoms  into  ridges  -  one  was 
parallel  to  the  x-axis  and  the  other  two  were  parallel  to  the 
y-axis.   See  Figures  10a  and  11a.   Despite  the  different  ways  the 
atoms  were  initially  piled  and  despite  the  fact  that  different 
random  number  seeds  were  used  for  the  two  ridges  paralleling  the 
y-axis,  there  was  a  remarkable  resemblance  in  the  gross  appearance 
of  equilibrium  microstates  of  each  experiment  -  these  are  sketched 
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as  Figures  10b,  lib  and  lie.   The  slots  appearing  in  these  sketches 
are  all  in  the  x  =  1  plane,  which  was  not  completely  filled  in  the 
initial  microstates,  i.e.,  the  holes  do  not  represent  several  atoms 
jumping  out  of  the  z  =  1  plane  but  the  rearrangement  of  a  hole 
originally  placed  there  after  it  has  received  some  additional  atoms 
from  the  higher  planes. 

Next,  four  simulations  were  run  in  which  the  atoms  were  piled 
into  valleys  of  different  step  widths.   The  initial  microstates 
are  sketched  in  the  "a"  parts  of  Figures  12  through  15.   Note  the 
similarity  between  Figures  10a  and  12a  and  between  Figures  11a  and 
l4a  -  in  effect,  since  periodic  boundary  conditions  were  applied, 
Figure  12a  may  be  obtained  from  Figure  10a  by  translating  Figure  10a 
one-half  the  active  volume  distance  in  the  y-direction,  and  simi- 
larly  for  Figures  11a  and  l4a  with  translation  in  the  x-direction. 
In  any  event,  it  may  be  expected  that  some  overall  gross  similarity 
of  the  final  microstates  may  appear.   Comparison  of  Figures  10b  and 
12b,  and  Figures  lib  and  l4b  tends  to  verify  this  expectation. 
Note  that  again  the  holes  appearing  in  the  final  microstates  were 
a  result  of  rearrangements  and  partial  fill-in  of  the  original 
holes  in  the  z    =   1   planes.   Figures  13a  and  l5a  indicate  valleys  of 
narrower  step  width,  and  final  microstates  achieved  are  sketched 
as  Figures  13b  and  15b. 

Since  there  was  a  preference  for  alignment  along  the  (lio)  di- 
rection, the  two  final  simulations  started  with  all  atoms  aligned 
in  that  direction.   One  was  visualized  as  a  truncated  (1  -1  1) 
plane  (Figure  l6a),  and  the  other  as  a  truncated  (-1  1  1)  plane 
(Figure  17a).   The  500th  microstates  are  sketched  as  Figures  l6b 
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and  17b.   Note  that  in  each  case  atoms  were  still  piled  high  and 
no  plane  was  completely  filled  -  these  facts  indicate  that  the 
system  probably  had  not  reached  equilibrium.   These  initial  micro- 
states  were  a  bit  more  unphysical  than  the  others  tested  since  they 
started  with  steps  several  atomic  layers  high  on  2  sides  and  one 
atomic  layer  high  on  the  other .   This  permitted  overhangs  when 
atoms  "jumped  off"  the  steep  sides,  but  this  was  not  considered 
seriously  detrimental  to  the  results. 
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IV.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.   CONCLUSIONS 

As    indicated   in    the    discussion    of   the  perfect    surface   simu- 
lation  results,    raost    "experiments"   yielded   results    which    imply 
that    the    (110)    orientation   was    most    natural    for    sketching    equilib- 
rium raicrostates.      Here,    "natural"   meant    that    there   were   fewer 
jagged  edges    and  protruding   atoms    than    if    the   microstate   were 
sketched   in    the    (100)    orientation.      Note    that    the    initial    con- 
figuration   sketches    (the    "a"  part    of   Figures    3    through    17)    are  all 
in    the    (100)    orientation,    since    this    was    how    the   crystal    was    initi- 
ally   visualized  and   developed.      The   following    discussion   will 
neglect    the    one   unphysical    case   achieved    -    that    of  an    initially 
perfect    surface   at    10,000   °K . 

The    first    important    result    of   testing   this    model   was    that    there 
was    one   thread   of   similarity   running    through   all    the   experiments    - 
in    the    final    microstates    achieved,    all    atoms    coalesced    into   steps 
either    parallel    or    perpendicular    to  a    (llO)    direction.      This   was 
most    encouraging   since    the    initial   state   had   no   effect    on    the   final 
step    direction.       It   was    further    encouraging   since    when   building   an 
fee,    (001)    crystal    from   marbles    and   glue,    atoms    on    the   surface   are 
more    closely    arranged   in   a    <110)    direction    than    in   a    <100)    direction, 
The   physical    reason    for    such   a  preference    lies    in    energy    considera- 
tions:   in   an    fee,     (001)    crystal,    NNl's    are    located  along    (110> 
directions,    and      NN2's    along    <100)    directions.      Therefore,    the 
lower    energy    state   sought    by   a    system    intuitively   and   deliberately 
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built  into  the  transition  probabilities  of  this  model  indicate  a 
preference  for  the  (llO)  direction.   Again,  it  is  a  pure  energy 
consideration  which  dictated  that  the  atoms  coalesce  and  not  dis- 
perse randomly.   The  minimum  energy  condition  was  achieved  for  as 
many  atoms  as  possible  staying  close  together.   Occasionally,  iso- 
lated surface  vacancies  or  stub  atoms  would  appear,  but  they  were 
quickly  repaired  or  moved  so  that  these  defects  tended  to  coalesce. 

The  second  conclusion  involved  the  final  shape  of  piled  atoms. 
Atom  piles  tended  to  be  stable  in  a  stepped  configuration  over 
several  microstates,  and  this  agreed  with  preconceived  notions  of 
surface  shape.   Figure  27  shows  a  sketch  due  to  D.E.  Harrison  [8j 
indicating  his  concept  of  a  surface  equilibrium  microstate.   Ac^ 
cording  to  the  model  developed  and  tested  here,  the  sketch  is 
essentially  correct  if  it  is  drawn  in  a  (110)  orientation  for  its 
stepped  features.   However,  this  model  indicated  that  isolated 
holes  and  stubs  do  not  persist  at  equilibrium.   Isolated  vacancies 
or  "holes"  were  as  mobile  as  isolated  stub  atoms.   Isolated  groups 
or  atoms  and  holes  appeared  to  move  slowly  as  a  group  to  coalesce 
with  larger  aggregates. 

The  third  observation  of  importance  involved  breaks  and  their 
repair  in  elongated  strings  of  surface  atoms.   Figure  28  shows  a 
series  of  sketches  illustrating  this  phenomenon.   When  the  motion 
of  the  atoms  led  to  the  movement  of  atoms  in  opposite  directions, 
a  local  "necking  down"  of  surface  aggregates  of  atoms  was  observed. 
Random  motions  led  to  selective  breaks  at  such  thin  steps,  but 
unless  the  atoms  in  the  two  resulting  steps  moved  several  atomic 
diameters  apart  the  break  actually  repaired  itself.   Conversely, 
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unless  the  break  repaired  itself  by  supplying  several  atomic 
diameters  of  step,  the  break  was  likely  to  recur  within  10  or 
20  microstates.   This  behavior  was  observed  repeatedly,  and  was 
not  really  surprising  although  it  was  not  expected  beforehand. 

B.   RECOMMENDATIONS 

First,  it  is  recommended  that  this  model  be  further  studied 
and  refined  to  obtain  a  more  realistic  picture  of  metal  surface 
microstates.   The  results  achieved  and  conclusions  drawn  here  must 
be  considered  only  as  preliminary  indications  since  this  was  a 
first  attempt  at  a  new  physical  model  and  since  the  testing  was 
not  extensive. 

Secondly,  it  is  recommended  that  prior  to  extensive  experi- 
mentation to  find  equilibrium  microstates,  the  problem  of  model 
dependent  transport  be  more  thoroughly  investigated.   Since  large 
aggregates  of  atoms  provided  the  most  interesting  results  here, 
this  testing  should  be  carried  out  with  groups  of  at  least  10 
atoms  on  a  perfect  surface.   It  is  envisioned  that  a  practical  test 
might  involve  calculating  the  center  of  mass  of  these  aggregates 
for  several  hundred  microstates  and  computing  the  average  motion 
of  their  center  of  mass.   If  trends  appear  (as  in  the  case  of  a 
single  stub  atom,  tested  here),  then  appropriate  corrections  must 
be  made  in  the  computer  program.   Perhaps  the  easiest  possible 
change  to  install  involves  the  sequence  in  which  the  atoms  are 
looked  at  as  a  microstate  is  generated.   As  it  was  used  here, 
every  microstate  was  generated  by  having  the  computer  investigate 
each  atom  in  order  by  lattice  site  number,  which  meant  that  rows 
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were  investigated  from  left  to  right  and  from  bottom  to  top  in 
each  z-plane.   It  may  be  more  realistic  to  generate  every  other 
raicrostate  in  a  different  way,  e.g.,  by  investigating  each  column 
from  bottom  to  top  and  from  left  to  right  in  each  z-plane. 

Third,  it  is  recommended  that  the  following  changes  be  made  in 
order  to  improve  the  model: 

a.  Rewrite  this  fee  program  to  generate  the  active  lattice 
volume  in  a  (110)  orientation.   This  would  require  a  new  lattice 
generator,  new  lattice  loading  steps,  new  NN  assignment  segments, 
and  recalculation   of  potential  energies. 

b.  Include  consideration  of  NN3's  and  NN^'s  -  this  would 
require  new  input  parameters  for  the  truncated  Morse  potential 
plus  NN3  and  NN4  assignment  segments. 

c.  Allow  floating  point  values  for  -flie  atom  positions  and 
introduce  kinetic  energy  for  the  atoms.   It  is  expected  that 
introducing  a  normal  distribution  of  kinetic  energies  would  re- 
duce the  total  energy  of  perfect  surface  atoms  sufficiently  to 
start  them  "jumping"  out  of  their  positions  without  resorting  to 
the  equivalent  temperature  idea  introduced  here.   This  would  al- 
low for  the  introduction  of  lattice  temperature  as  an  input  para- 
meter and  controlling  factor  in  lattice  activity. 

d.  Write  a  lattice  loader  that  would  generate  a  more  realistic 
lattice  -  e.g. ,  one  in  which  the  atoms  were  randomly  placed  sub- 
ject to  physical  restraints  such  as  having  at  least  one  NNl  and 
occupying  nearly  perfect  lattice  positions  (i.e.,  this  assumes 

that  floating  point  lattice  positions  have  already  been 
incorporated) . 
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e.      Allow    for    vacancies,     inter st it ials   and    impurities.      This 
would   require   consideration    of   lattice   relaxation   and   strain 
effects,    and  would   require   new  potentials    for    impurity-copper    and 
impurity-impurity    interactions. 
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APPENDIX  A 
THE  DETAILS  OF  THE  SIMULATION  MODEL 

A.   ASSUMPTIONS 

Perhaps  the  greatest  assumption  made  in  this  calculation  is 
that  metal  crystal  binding  energies  can  be  computed  from  a  single 
composite  potential  energy  function,  with  the  interactions  between 
atoms  considered  as  independent  two-body  problems.   The  essential 
feature  of  this  assumption  was  to  treat  the  atoms  of  a  crystal  lat- 
tice as  quasi-hard  spheres  which  interact  by  some  ^'unknown"  mecha- 
nisms to  yield  observed  crystal  behavior.   The  "unknown"  mechanisms 
include,  in  the  general  case,  attraction  due  to  ionic,  covalent, 
metallic  and  van  der  Waals  forces,  and  repulsion  due  primarily  to 
electron  cloud  overlap.   The  specific  mechanisms,  however,  were  of 
no  concern.   No  quantum  mechanical  considerations  were  attempted  - 
that  is,  the  calculations  involved  were  all  classical  once  the  ap- 
propriate composite  potential  function  was  assumed.   First  and 
second  nearest  neighbor  interactions  were  considered. 

The  only  metal  visualized  was  copper ,  which  forms  a  face- 
centered  cubic  (fee)  lattice.   For  simplicity,  the  (001)  crystal 
orientation  was  assumed.   Lattice  sites  were  visualized  as  ordered 
real  triples  in  the  positive  octant  of  a  rectangular  coordinate 
system.   The  coordinate  planes  and  axes  bordering  the  positive 
octant  were  considered  as  belonging  to  the  positive  octant.   See 
Figure  1.   Although  a  rectangular  parallelepiped  was  always  gene- 
rated to  provide  possible  lattice  locations  (hereafter  called  the 
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"active  lattice  volume"),  appropriate  changes  in  the  lattice  gene- 
rator insured  that  the  initial  configuration  of  the  surface  could 
take  any  desired  shape  with  the  sole  precondition  that  the  only 
sites  available  were  perfect  lattice  sites.   The  active  lattice 
volume  was  considered  as  merely  a  framework  of  sites  that  were 
available  for  atoms  -  whether  the  sites  were  actually  occupied  or 
not  depended  on  the  details  of  the  micr estate. 

It  was  further  assumed  that  atomic  vibrations  about  lattice 
sites  could  be  neglected,  so  there  were  no  zero-point  energy  con- 
siderations or  fluctuations  in  lattice  potential  due  to  atomic 
motion.   Lattice  relaxation  due  to  crystal  vacancies  was  ignored, 
as  were  interstitial  atoms  and  their  consequent  lattice  strain 
effects.   Furthermore,  no  impurity  atoms  were  allowed.   With  these 
assumptions,  it  was  possible  to  describe  each  available  lattice 
site  as  ordered  integer  triples  and  use  a  single  composite  po- 
tential energy  function  for  every  two-body  pair  in  the  system. 
This  simplified  the  model  considerably  and  yielded  a  great  re- 
duction in  the  potential  energy  calculations  required. 

In  addition,  it  was  assumed  that  only  the  top  several  atomic 
layers  of  the  surface  would  provide  interesting  results  -  in 
effect,  this  meant  that  the  active  lattice  sites  generated  were 
placed  on  top  of  a  perfect  lattice  substructure  which  was  not 
allowed  to  interact  with  the  surface  sites  in  any  way  other  than 
to  be  used  in  calculating  surface  site  potential  energy;  this  fur- 
ther implied  that  above  the  active  lattice  sites  all  locations  were 
always  kept  vacant.   Specifically,  lattice  sites  described  by  2<0 
were  considered  always  full  -  see  Figure  1.   Lattice  sites 
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described   by   0  ^    2   ^    z  were    considered   as    available  for    con- 

'  max 

taining  an  atom,  depending  on  the  particular  microstate  initialized 

or  achieved  during  a  calculation.   Lattice  sites  described  by 

2.  ">   2.  were  considered  always  empty. 

max  ^ 

Finally,  this  model  was  restricted  to  even  numbers  of  planes 
in  each  of  the  principal  coordinate  directions  to  simplify  the 
calculations  involved  in  assigning  first  and  second  nearest  neigh- 
bors to  each  available  site.   This  means  that  x    ,  y    ,  and  z 

max    max        max 

must  all  be  odd,  since  each  coordinate  plane  was  included  as  part 
of  the  active  lattice. 

B.   THE  COMPUTER  PROGRAM 

1.   The  Lattice  Generator 


a).  Logic.   Before  any  calculations  could  be  made,  the  com- 
puter had  to  be  told  what  the  available  sites  for  atoms  were  and 
whether  or  not  these  sites  had  an  atom  -  that  is,  whether  a  given 
site  was  occupied  or  vacant.   For  the  (001)  orientation  of  an  fee 
crystal  (hereafter,  "fee,  (001)  crystal"),  the  essential  technique 
involved  relied  on  the  fact  that  the  sum  of  the  site  coordinates 
must  be  an  even  (or  an  odd)  integer  (this  point  is  particular  to 
this  lattice  shape  and  orientation).   Refer  to  Figure  29.   For 
example,  the  Cartesian  point  (1,  1,  1)  has  a  sum  of  coordinates 
equal  to  three,  an  odd  integer  -  therefore,  this  point  is  not  an 
acceptable  lattice  location.   On  the  other  hand,  Cartesian  points 
(1,  0,  1)  and  (1,  1,  0)  each  have  a  sum  of  coordinates  equal  to 
two,  an  even  integer  -  these  points  are  acceptable  lattice  locations 
and  are  assigned  integer  labels  when  they  are  reached  in  the  lattice 
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generator  segment  of  the  program.   Figure  30  shows  the  numbering 
scheme  for  lattice  points. 

The  lattice  was  generated  by  a  set  of  nested  DO  loops. 
One  plane  at  a  time  was  generated,  and  these  planes  were  parallel 
to  the  x-y  plane.   The  coordinates  of  each  acceptable  site  and  its 
occupation  index  were  stored  as  vector  arrays.   An  "occupation 
index"  (vector  "NOCC(I)")  merely  indicated  whether  or  not  a  site 
was  occupied:  an  index  of  zero  implied  that  the  site  was  vacant, 
while  an  index  of  one  meant  that  it  was  occupied.   No  values  other 
than  zero  or  one  were  allowed  for  occupation  indices. 

b) .  The  Initial  Microstate.   In  order  to  provide  adequate 
testing  of  the  simulation  model,  the  initial  microstate  must  be 
easily  changed.   This  was  achieved  by  substituting  program  segments 
containing  one  or  more  cards  into  the  program.   Initial  microstates 
actually  used  were: 

1)  Perfect  surface 

2)  Perfect  surface  plus  extra  atoms 

(a)  randomly  dispersed  extra  atoms  -  Figures 
3a  and  4a 

(b)  a  monatomic  step  -  Figure  5^ 

(c)  pyramids  -  Figures  6a  and  7a 

(d)  pyramid  holes  -  Figures  8a  and  9a 

(e)  ridges  -  Figures  10a  and  11a 

(f)  valleys  -  Figures  12a,  13a,  l4a  and  15a 

(g)  (1  -1  1)  plane  -  Figure  l6a 
(h)   (-111)  plane  -  Figure  17a 

The  actual  program  segments  which  generate  these  different  initial 
conditions  are  listed  in  Appendix  B. 
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A  site  labeled  ^'IFUL"  was  used  to  represent  all  sites 

in  the  crystal  substrate,  that  is,  all  sites  for  z  <  0;  since  these 

sites  were  always  full  (occupied),  NC)CC(IFUL)  was  always  one.   A 

site  labeled  "IVAC"  was  used  to  represent  all  sites  above  the 

active  lattice  volume,  that  is,  for  2  >  z    ;  since  these  sites 
'         '  max 

were  always  vacant,  NOCC(IVAC)  was  always  zero.   These  assignments 
were  made  immediately  after  the  lattice  generator  and  were  based 
on  the  last  value  of  "M"  reached  in  the  lattice  generator. 
2.   Periodic  Boundary  Conditions 

a).  The  need.   Computer  storage  and  time  limitations  made 
impractical  the  use  of  an  arbitrarily  large  active  lattice  volume. 
A  practical  lattice  size  involved  at  most  a  few  thousand  active 
lattice  sites.   This  required  that  the  active  volume  size  be  limit- 
ed to  a  Cartesian  space  of  roughly  20  x  20  x  10,  which  corresponds 
to  2,000  active  lattice  sites  and  a  copper  crystal  measuring  about 
50  A  X  50  A  X  25  A.   Such  a  small  crystal  with  six  surfaces  com- 
plicated the  study  of  surface  effects. 

b) .  The  method.   One  surface  (the  bottom  of  the  active 
lattice  volume)  was  eliminated  by  assuming  that  the  active  lattice 
volume  was  superimposed  on  a  perfect  crystal  substrate  which  was 
never  allowed  to  vary.   Four  other  surfaces  (the  sides  of  the 
active  volume)  were  eliminated  by  assuming  periodic  boundary  con- 
ditions which  resulted  in  an  essentially  infinite  single  surface 
(the  ^'top"  of  the  active  volume)  that  varied  periodically  in  con- 
figurationo   This  infinite  surface  may  be  visualized  as  part  of  a 
crystalline  surface  far  from  a  grain  boundary. 
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The  overall  scheme  was   represented  as  follows  (see 
Figure  31):  Each  active  volume  plane  was  laterally  surrounded  by 
exactly  similar  active  volume  planes  that  had  been  translated  in 
the  X  and/or  y  dir ection( s ) .   In  Figure  31a,  the  actual  active 
volume  was  visualized  as  the  center  square;  the  remaining  eight 
squares  were  the  "boundary  volume  planes"  which  had  exactly  the 
same  configuration  as  the  center  square.   Square  2  was  achieved 
by  translating  square  1  a  full  microcrystal  distance  in  the  "-x" 
direction,  and  square  3  was  achieved  by  translating  square  1  a 
full  microcrystal  distance  in  the  "+x"  direction.   Squares  4-5-6 
were  achieved  by  the  similar  process  of  translating  squares 
2-1-3  a  full  microcrystal  distance  in  the  "+y"  direction,  and 
similarly  for  squares  7-8-9. 

The  active  volume  z-plane  of  Figure  31b  is  por- 
trayed as  having  six  x-planes  and  six  y-planes,  and  a  total  of  l8 
lattice  sites;  three  of  these  lattice  sites  are  labeled  (A,B,C), 
and  the  location  of  the  remaining  sites  are  indicated  by  X's.   The 
images  of  the  labeled  sites  are  indicated  in  each  of  the  surrounding 
boundary  planes.   Without  the  boundary  planes,  a  vector  [l,l,o] 
from  active  volume  site  "C"  would  yield  a  point  in  free  space. 
However,  with  the  boundary  planes,  as  in  Figure  31b,  that  same 
vector  from  "C"  yields  a  site  labeled  "B"  in  the  boundary  plane. 
In  effect,  looking  to  the  right  from  the  active  volume  plane,  one 
passes  through  infinity  and  comes  upon  the  left  edge  of  the  same 
active  volume  plane  in  the  same  row  from  which  one  started.   Simi- 
larly, a  vector  [o,6,oJ  from  site  "A"  brings  one  right  back  to 
site  "A".   (This  is  merely  an  illustration,  and  actual  crystals 
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generated  typically  had  20  or  so  planes  in  the  x  and  y  directions, 
so  that  returning  to  the  same  site  was  not  so  easy  in  practice.) 
During  the  calculation,  any  image  sites  reached  were  treated  as 
active  volume  sites.   Specifically,  if  in  the  active  volume,  sites 
A  and  C  were  occupied,  but  site  B  was  vacant,  then  in  all  eight 
boundary  planes,  sites  A  and  C  were  occupied  and  site  B  was  vacant. 
This  gave  a  periodic  character  for  site  occupation  as  one  moved 
across  a  given  z-plane,  hence  the  name  periodic  boundary  conditions. 
In  actuality  only  one  active  volume  was  used,  but  the  ""image  site" 
and  "boundary  plane"  ideas  explained  here  help  to  illustrate  how 
periodicity  was  achieved. 

c) .  Second  nearest  neighbor  assignment  segment.   The  NN2 
assignment  segment  is  easier  to  understand  than  the  NNl  assign- 
ment segment,  and  so  is  considered  first.   Figure  19  illustrates 
the  positions  of  the  six  NN2's  in  a  fee,  (001)  crystal,  using  the 
numbering  scheme  developed  for  the  program.   To  find  the  site 
number  of  any  second  nearest  neighbor  of  a  given  atom,  the  given 
atom  was  visualized  as  the  center  of  coordinates  of  a  small  mobile 
coordinate  system  as  depicted  in  Figure  19.   In  general,  the  site 
at  the  origin  was  labeled  "I",  and  the  six  positions  for  NN2 ' s  were 
labeled  "J",  where  J  varied  from  one  to  six.   A  double  subscript 
notation  was  chosen  and  named  array  "NBRTWO( J , I) "  -  a  mnemonic  for 
"neighbor  two"  -  where  subscript  J  refers  to  one  of  the  six  NN2 
positions,  subscript  I  refers  to  the  label  of  the  lattice  site 
being  investigated  (given  to  that  site  by  the  lattice  generator), 
and  the  value  stored  at  the  particular  position  was  the  label  of 
the  lattice  site  at  position  J  with  respect  to  site  I. 
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In    the   general    case,    NfN2   position    one   was    labeled 
site   number    (I-l),    and  NN2   position    two   was    labeled   site   number 
(I+l).      NN2  positions    three   and  four    required   that    a    shift    of    two 
rows    of   atoms    in    the   "-y"   and    "+y"   directions   be   made.      For    even 
numbers    of   x  planes,    that    is,    for    IX   even,    the   appropriate   number 
of  site   numbers    to    "shift"    is    IX,    since   each   row   of    sites    in   a 
given  plane   contains    IX/2    sites    (IX   even).      Consequently,    NN2 
position    three   was    labeled   site   number     (I-IX)    and  NN2   position    four 
was    labeled  site    number    (I+IX).      NN2   positions    five   and   six   re- 
quired a    shift    of    two  planes    of   atoms    in    the   "-z"   and   "+2"   di- 
rections.     Since   each   row    of   a    given   plane   contained   IX/2    sites 
(IX    even),    and   there  were   lY   such   rows   per    plane,    each  plane   con- 
tained   (IX/2)*(IY)    sites.      Therefore,    two   such  planes    contain 
IX*IY   sites.      As    a    result,    NN2   position    five   was    labeled   site 
number    ( I-IX*IY)    and  NN2   position    six   was    labeled   site   number 
(I+IX*IY) . 

For    sites    on,    or    one   atomic    layer    in   from,    the   active 

volume   surface,    proper    labeling    of   lSfN2    sites    helped   achieve    the 

desired   effect    for    a    single    infinite   surface    varying   periodically. 

For    example,    a   site   with   NX(I)=0    (on    the    left    face   of    the  active 

volume)    or    NX(I)=1    (one   atomic    layer    in   from    the   left    face    of   the 

active   volume)    looking    for    its    neighbor    at    NN2   position    one   must 

"see"   that    neighbor    (call    it    I  '  )    with  NX(lM=2         -1    or    NX(I')=z         , 
^  '  ^       '       max  ^       '      max 

and  NY(I' )=NY(I) ,    and   NZ( I ' )=NZ( I ) .       In    the   site   labeling   scheme 
chosen,    this    represents    a    shift    of    (IX/2)-l    site   numbers    (IX   even). 
Similar    shifts    must    be   made   when    "looking    through    infinity"   from 
the   other    three   sides    of    the   active    lattice    volume,    and   this    was 
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essentially  how  the  imposed  periodic  boundary  conditions  were 
achieved  in  practice. 

For  sites  with  NZ(I)=0  or  NZ(I)=1  looking  for  their 
neighbor  at  NN2  position  five,  site  IFUL  (the  always  full  sub- 
strate position, was  reached.   Similarly,  for  sites  with  NZ(I)=2 

^  '  -'  '  ^     '      max 

or    NZ{I)=z         -1    looking    for    their    neighbor    at    NN2   position    six, 
^     '      max  ^  r-  7 

site    IVAC    (the   always    empty   above    lattice  position)    was    reached. 
In    this    way   a   perfect    substrate   below   the   active    volume   and   a    void 
above    the  active    lattice   were   maintained.      Further    specific    de- 
tails   of -frie   NN2    assignment    segment    are    omitted   here,    but    may   be 
found  by    studying    the   actual   program   between   statements    number 
195  and  201. 

d) .    First    nearest    neighbor    assignment    segment.      The   NNl 
assignment    segment    was   approached   in   much   the   same   way   as    the   NN2 
assignment    segment.      The  problem  was    complicated   by    the    fact    that 
there   are   12   NNl's    for    an    fee,    (001)    crystal,    and   there  are   22 
special    cases    (6    faces,    12    edges,    and   4    corners)    to   consider,    as 
opposed   to   six   NN2's    and   six   special    cases    for    the   NN2   assignment 
segment.      Figure    l8    illustrates    the  positions    of   the    12   NNl's, 
using    the   numbering    scheme    developed   for    the  prograiil.      Again,    a 
double    subscripted  array   was    chosen    to   store  NNl    information,    and 
it   was    called   "NBRONE( J, I) "    -   a   mnemonic    for    "neighbor    one"    -   where 
subscript   J    refers    to   one    of    the   12   NNl   positions,    subscript    I   re- 
fers   to   the    lattice   site   under    consideration,    and   the   value   at    the 
particular    array   position    is    the   site    number    of    the   site   at   position 
J   with   respect    to   site   I. 
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The  initial  complication  that  arose  was  an  "even-odd" 
dependence  on  the  x-value  of  the  site  for  8  out  of  12  NNl  positions 
Taking  this  into  account,  the  overall  problem  was  solved  by  first 
assigning  each  site  its  12  NNl's  without  regard  to  the  22  special 
cases.   These  were  computed  using  a  shift  technique  similar  to 
that  developed  in  the  NN2  assignment  segment.   The  22  special  cases 
involved  sites  on  the  faces,  edges,  and  corners  of  the  active  lat- 
tice volume,  and  required  particular  attention  in  order  to  appro- 
priately achieve  the  desired  periodic  boundary  conditions.   The 
method  used  then  involved  testing  for  each  of  the  22  special  cases. 
As  each  special  case  was  found  a  series  of  from  four  to  nine  ad- 
justments were  made  from  the  "general  NNl  sites"  already  stored. 
Sites  IFUL  and  IVAC  were  employed  in  a  manner  similar  to  that 
involved  in  the  NN2  assignment  segment.   Further  details  of  the 
calculation  are  tedious  and  will  not  be  examined  here.   The  actual 
procedure  may  be  gleaned  from  the  program  itself,  beginning  after 
the  lattice  generator  and  running  through  statement  195- 

Thus,  the  periodic  boundary  conditions  imposed  on  the 
active  lattice  volume  were  achieved  in  practice  via  the  use  of  the 
NNl  and  NN2  assignment  segments.   These  segments,  as  written,  re- 
quired that  IX,  lY  and  IZ  all  be  even.   This  was  not  considered  a 
very  severe  restriction,  and  there  is  no  a  prior  reason  to  expect 
that  the  results  were  biased  due  to  a  choice  of  even  numbers  of 
planes  in  the  coordinate  directions. 
3.   Potential  Energy 

a).  The  Composite  Potential  Energy  Function.   The  model 
chosen  assumed  that  the  crystal  may  be  approximated  by  a 
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latticework   of   quasi-hard  spheres,    so  a    scheme    for    holding    indi- 
vidual   atoms    together    was    needed.      The  problem    is    essentially   an 
N-body   problem,    and    therefore    no   general    exact    solution    exists. 
However,    approximating    the  potentials    by    an    independent    set    of 
two-body   potentials    yields   a    solvable   approximation    to    the  physical 
system,    and  was   a    standard   technique   used    in    radiation    damage 
studies . 

Since   the   approach    taken    involved    the  solution    of 
independent    two-body  potentials,    the   next    task    involved  a    judicious 
choice  of  potential    function.      The   choice   began  with    the  Morse 
Potential,    proposed   in    1929   as    an   appropriate  potential    for    the 
wave   equation   of   a    diatomic   molecule   LI3]: 

<^"^^        ^Morse"   ^  exp( -2a (r -r ^) )    -   2D  exp ( -a(r -r ^) ) , 

where 

D   =   dissociation  energy    of   the   molecular    bond 

r      =   equilibrium    interatomic    separation 

r    =    interatomic    separation 

a    =   constant   related   to   the   curvature  about    the 

potential    minimum   of    the   function 

Note   that    the  first    term   dominates    for    small    "r",    and   that    the 
second  term   dominates    for    large    "r". 

Girifalco   and  Weizer    calculated  Morse  potential    con- 
stants   for    15   cubic   metals,    and    these  have   found   satisfactory 
agreement    with   experiment   C?] •      Anderman    used   essentially    the   same 
procedure   as   Girifalco  and   Weizer,    but    used  a    truncated  potential 
and   obtained  parameters    for    a    copper    Morse-type  potential    for    two, 
three,    and   four    nearest    neighbors    Cl] .      This    "truncated  Morse" 
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potential    with  parameters    for    two   nearest    neighbors    was    used  for 

the  range   r      <   r   <    r    ,    where   r      is   an    inner    cutoff   value   chosen 
o  C  B 

as    a    lower    bound   for    the   truncated  Morse  potential,    and  r       is    the 
truncation    distance. 

GGMV  arrived   at   reasonable   short    range   results    using   a 
Born-Mayer    type   ejqjonential   approximation    for    the    interatomic 
potential : 

(A-2)      ^GibsonII=    D^-P(-M^-^o))' 

where  D,  a,  r  and  r   have  the  same  meanings  as  in  V  . 

'   '        o  ^         Morse 

This  "Gibson  II"  potential  was  found  to  be  a  good  approximation 

for  small  interatomic  separations,  and  so  was  used  for  short 

distances,  0  <  r  <  r.,  where  r   is  an  outer  cutoff  value  chosen  as 
A         A 

an  upper  limit  for  the  region  of  Gibson  II  potential  validity.   In 

general,  r  f   x      because  these  two  potentials  do  not  match  smoothly 
A     B 

-  see  Figure  20.   The  "mismatch  region"  was  fitted  by  a  series  ex- 
pansion with  a  cubic  function,  since  only  four  pieces  of  information 
were  available  to  match  as  boundary  conditions  (these  were  the  value 
and  slope  of  the  Gibson  II  potential  at  r  ,  and  the  value  and  slope 

of  the  truncated  Morse  potential  at  r  )  -  the  most  general  cubic 

B 

(A-3)   Ar^  +  Br   +  Or  +  D 

has  four  unknowns,  which  just  equals  the  information  available. 
The  cubic  is  a  standard  choice  for  this  "linking  potential". 

Solving  for  the  coefficients  of  the  cubic  potential 
involved  solving  four  simultaneous  linear  equations  in  four  un- 
knowns.  The  technique  of  solution  involved  here  required  putting 
the  coefficients  of  the  four  linear  equations  into  a  matrix  which 
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was  then  transferred  to  Subroutine  CROSYM  (developed  at  Lockheed) 

which  performed  manipulations  on  the  matrix  (using  the  "Method  of 

Crout").   An  answer  vector  was  coupled  back  to  the  MAIN  program. 

The  actual  potential  function  used  in  the  program  was 

a  composite  potential,  employing: 

the  Gibson  II  potential,  for  0  <  r  <  r  , 

the  cubic  potential,  for  r  ^  r "^  r  , 

A  n 

and  the  truncated  Morse  potential,  for  r   <  r  <  r  . 
Figure  21  shows  the  overall  composite  potential.   Note  that  all 
interatomic  separations  used  in  this  program  were  greater  than  r  , 
that  is,  the  cubic  and  Gibson  II  potentials  were  never  used,  in 
practice.   However,  they  were  included  to  enchance  program  flexi- 
bility, since  it  was  anticipated  that  future  improvements  of  this 
model  would  allow  atomic  vibrations  which  may  require  the  use  of 
these  other  potentials.   All  these  formulas  are  for  interatomic 
spacings  in  lattice  units.   A  lattice  unit  is  the  distance  between 
two  NNl's  in  a  single  coordinate  direction. 

b) .  Site  Potential  Energy.   The  next  step  was  to  assume 
that  each  lattice  site  experienced  this  composite  potential  as  the 
appropriate  two-body  interatomic  interaction  potential.   Program 
statements  1050  through  1100  set  the  square  of  the  equilibrium 
interatomic  separation  with  the  simple  formula 

(A-4)   NNDIS2(N)  =  2*N. 
This  formula  is  exact  for  the  fee,  (001)  crystal  for  NNl  through 
NN13  in  lattice  units.   Using  these  values  for  distance  squared,  a 
potential  was  computed  for  each  ith  nearest  neighbor  (i=l,2)  based 
on  its  distance  from  any  particular  site.   Consider  the  sketch  in 
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Figure  21.   The  program  computed  and  stored  V(r  )  and  V(r  )  .   It 
then  looked  at  each  site  and  (whether  it  was  occupied  or  not),  it 
computed  an  energy  associated  with  that  site  based  on  nearby  oc- 
cupied sites  out  to  and  including  the  second  nearest  neighbor  sites 
This  energy  was  then  multiplied  by  PFIV  (a  mnemonic  for  "point 
five"  and  equal  to  0.5  typically,  but  included  as  a  convenient 
parameter  for  calculation)  since  on  the  average  about  half  the 
potential  was  associated  with  each  atom  of  the  many  two-body  pairs 
envisioned.   Values  of  PFIV  other  than  0.5  gave  a  different  value 
of  effective  lattice  temperature,  and  in  practice  lattice  tempera- 
ture was  varied  by  changing  PFIV.   The  true  equivalent  lattice 
temperature  based  on  potential  energy  only  was  defined  by 

(A-5)    EQTEMP  =  0.5*TEMP/PFIV, 
where  EQTEMP  is  the  equivalent  lattice  temperature 

TEMP  is  the  input  temperature  used  (=1000.0  °K) . 
Figure  22  is  a  log-log  plot  of  (A-5) • 

In  this  way  each  site  was  given  a  potential  energy 
based  on  nearby  occupied  sites.   These  energies  were  stored  as  the 
vector  PE(I),  where  I  was  the  site  number.   An  energy  associated 
with  both  occupied  and  vacant  sites  was  required  in  order  to  deter- 
mine relative  transition  probabilities  from  occupied  lattice  sites 
to  nearby  vacant  sites  -  and  this  transition  probability  was  re- 
lated to  the  energy  change  between  sites  by  what  was  essentially  a 
Boltzmann  factor . 

4.   Transition  Probability 

a).  Definition.   If  a  given  lattice  site  was  occupied  and 
had  one  or  more  vacant  sites  as  NNl's,  it  had  (statistically)  a 
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finite  probability  of  moving  from  its  present  site  to  the  nearby 
vacant  site.   In  practice,  when  an  occupied  lattice  site  was 
reached,  an  energy  difference  array  was  defined  between  the  site 
reached  and  each  of  its  12  NNl's.   The  array  was  defined  as 
"DE(J,I)"  -  a  mnemonic  for  "difference  in  energy"  -  where  I  was  the 
occupied  site  reached  and  J  (running  from  1  through  12)  referred  to 
the  position  of  one  of  the  NNl's  for  site  I.   The  value  stored  at 
the  array  position  was  the  energy  difference  between  site  I  and  the 
Jth  NNl  position  with  respect  to  site  I;  that  is,  if  K  were  the 
site  number  for  the  Jth  NNl  position  with  respect  to  site  I,  then 

(A-6)      DE(J,I)  =  PE(K)  -  PE(I). 
Since  the  crystal  was  always  bound,  every  PE(I)  and  PE(K)  was 
negative,  which  implied  net  attraction  or  crystal  cohesiveness . 
However,  i:£(J,I)  may  have  either  positive  or  negative  values,  de- 
pending on  whether  site  I  or  site  K  was  more  tightly  bound.   A 
negative  DE(J,I)  implied  that  site  K  was  more  tightly  bound  than 
site  I,  while  a  positive  DE(J,I)  meant  that  site  I  was  more 
tightly  bound  than  site  K. 

Complications  arose  since  it  was  unphysical  to  expect 
that  any  atom  (occupied  site)  had  no  NNl's.   No  atom  was  permitted 
to  jump  into  a  site  that  was  already  occupied  -  this  implied  destruc- 
tion or  loss  of  matter,  which  was  not  allowed.   The  case  of  inter- 
changing the  atoms  between  two  adjacent  occupied  sites  was 
physically  allowable,  but  uninteresting  since  it  did  not  change  the 
shape  of  the  microstate,  and  so  was  not  considered.   Finally,  tran- 
sition from  a  vacant  site  was  not  permitted,  since  this  implied 
creation  or  addition  of  matter  which  was  not  considered. 
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sites    with   no   NNl's    were   not    encountered    in   practice 
(exept    for    the   unphysical    "experiment"  at    10,000      K) .       The   problem 
of    jumping    from   a    vacant    site   and    "creating"   an   atom  was    easily 
solved  by    skipping   further    consideration    of   a    site    if    it    were 
vacant.      The    cases    of   jumping    from   one   occupied   site    to   another 
(matter    destruction)    and   of    interchanging    the   atoms    in   adjacent 
occupied   sites    (uninteresting    interchange)    were    taken    care    of   by 
defining   an   arbitrarily    large   DE(J,I)    =    1000.0    in    the  program 
whenever    an    occupied   site    (I)    found   that    the   site   at    its    Jth  NNl 
position   was    occupied.      Later ,    when    transition   probabilities    were 
computed,    any    DE(J,I)    ^    900.0  was    ignored,    that    is,    no   transitions 
were  allowed    (or    in    other    words,    the    transition   was    given   a  proba- 
bility   of   occurrence   of   0.0).      Sites    in   the    layer    z    =   z 

'  ^         max 

"looking"  out  the  top  of  the  active  lattice  volume  always  "saw" 
site  IVAC  which  was  not  allowed  to  be  occupied  -  such  sites  were 
also  given  a  DE(J,I)  =  1000.0  to  avoid  atom  loss. 

Intuitively,  physical  systems  reside  at  or  near  the 
bottom  of  potential  wells,  so  one  expects  that  for 

DE(J,I)  <  0 

the  atom  at  site  I  would  probably  jump  to  site  K,  since  it  would 
then  have  a  lower  potential.   In  fact,  by  the  way  transition  proba- 
bility was  defined,  this  was  exactly  what  happened.   The  transition 
probability,  p,  was  defined  as  [lO]: 

(A-7a)     p  =  1  DE(J,I)  ^  0 

(A-7b)  p    =    exp(-DE(J,I)/kT)  DE(J,I)    >    0 

where  k   =    Boltzmann's    constant 

T    =   absolute    temperature. 
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Note  that  in  the  program,  "kT"  was  defined  at  "TAU" ,  that  is, 
(A-8)     TAU  =  k*T. 

For  two  or  more  transitions  which  satisfy  condition 
(A-7a),  the  computer  searched  for  the  deepest  well,  and  carried 
out  the  transition.   Here,  that  would  mean  that  for  DE(J,I)  ^  0 
that  the  atom  at  site  I  jumped  to  site  K.   "Carry  out  the  transi- 
tion" meant  that  the  computer  zeroed  the  occupation  index  of  site 
I  and  placed  a  "1"  in  the  occupation  index  of  site  K.   In  the 
event  that  two  or  more  wells  had  the  same  (or  zero)  depth  for  the 
deepest  nearby  well,  the  computer  chose  one  of  the  sites  by  a 
random  number  process.   This  process  involved  matrix  FTEST  which 
was  loaded  by  nesting  DO  loops  at  the  beginning  of  the  program. 
Matrix  FTEST(I,J)  was  12  x  12  and  had  elements 

0.0  below  the  diagonal, 

1.0  on  the  diagonal, 

FLOAT(I)/FLOAT(J)     above  the  diagonal. 
When  the  energy  difference  was  the  same  for  a  jump  to  two  or  more 
different  sites,  the  number  of  sites  with  the  same  energy  was  used 
as  J  and  a  DO  loop  index  used  to  vary  I  from  1  to  J.   A  random 
number  was  chosen  and  the  FTEST(I,J)  values  compared  with  the 
random  number.   If  the  FTEST(I,J)  value  equalled  or  exceeded  the 
random  number,  the  transition  was  made  to  the  site  which  was  the 
Ith  one  which  had  that  particular  negative  energy  difference.   If 
the  random  number  exceeded  the  FTEST(I,J)  value,  the  loop  index 
was  incremented  and  the  new  FTEST(I,J)  value  compared  with  the 
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same  random  number,  and  so  on,  until  one  FTEST(I,J)  value  equalled 
or  exceeded  the  random  number.   In  any  case,  a  transition  was 
always  made  for  negative  (or  zero)  energy  differences. 

The  computer  decided  whether  or  not  two  energy  dif- 
ferences were  "the  same"  by  using  a  small  parameter   "EPSLON"  -  a 
mnemonic  for  "^epsilon".   Two  energies  were  considered  to  be  "the 
same"  if  they  were  within  EPSLON  of  each  other.   In  practice, 
EPSLON  =  1.0  X  lO"^. 

For  situations  in  which  all  DE(J,I)  were  positive,  the 
transition  decisions  were  based  on  a  random  number  process.   First, 
the  smallest  positive  energy  difference  was  found.   Then  a  proba- 
bility for  this  lowest  positive  energy  difference  was  computer  by 
formula  (A-7b) .   Next  a  random  number  was  chosen  on  the  interval 
(0,1).   If  the  transition  probability  equalled  or  exceeded  the 
random  number,  the  computer  carried  out  the  transition,  and  moved 
on  to  the  next  occupied  site  to  repeat  the  whole  process.   If  the 
probability  did  not  equal  or  exceed  the  random  number,  the  com- 
puter found  the  next  lower  positive  energy  difference,  calculated 
a  probability,  and  compared  it  to  a  new  random  number.   This  process 
was  continued  until  either  a  transition  was  made  to  a  state  of 
higher  energy  or  the  site  under  study  ran  out  of  first  nearest 
neighbor  vacant  sites  -  in  this  latter  case,  no  transitions  were 
made  to  any  other  site,  and  the  computer  moved  on  to  the  next 
occupied  site  to  repeat  the  procedure  for  the  new  site.   During 
the  calculation,  if  two  or  more  sites  had  the  same  positive  energy 
difference,  these  sites  were  treated  sequentially  in  a  manner  simi- 
lar to  sites  of  the  same  negative  energy  differences.   Positive 
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site   energies    were  considered   "the   same"   if   they   were  within 
EPSLON   of   each    other    -    this    was    essentially    the   same  procedure 
used   for    the   negative   DE(J,I)'s. 

After    any    transition,    the   potentials    of   the   siteSv 
affected  were    "patched".      The   neighbors    of   the    "old   site",    that 
is,    the  previously    occupied   one,    lost    an   NNl    or    an   ^fN2 ,    therefore, 
PFIV   times    the  proper    energy   was    subtracted  from   each    of    these 
NNl's    and  NN2's.       In  addition,    the   neighbors    of   the    "new   site" 
each   gained  an    NNl    or    an    NN2 ,    so  PFIV    times    the   proper    energy   was 
added   to   each   of    these   NNl's    and  NN2's. 

The    question    of    how  appropriate  are    these  probabilities, 
(A-7a)    and    (A-7b) ,    may   arise.      Such  probabilities    have   been    used 
in    the    literature    on   crystal    growth,    and  have   been   shown   by    sta- 
tistical  arguments    to  be   correct    for    an    infinitely    long   chain    of 
microstates.      According    to   Leamy   and  Jackson: 

"averages    over    a    chain   of   microstates    obtained 
Cthrough    the   use    of   these    transition   probabilities] 
...    will    converge    to    the   equilibrium  properties    of 
the   system    in   the    limit    of    infinite   chain    length. "[lo] 

b) .      Pseudo-Random  Number    Generator.      The    technique   used 
to  generate   random   numbers    was    a    standard   one    called    "pseudo- 
random  number    generation   by    the  multiplicative   congruential    method". 
The   term    "psuedo-random"    implied  that    the   numbers    generated  were 
not   really   random,    and   this    was,    in    fact,    true.      However,    such 
methocfeare   generally   accepted  provided   they  pass    certain   statistical 
tests    which   determine   whether    the   numbers    generated  are,  ±>x    all 
practical   purposes,    random.      The    tests    require   that    the   numbers 
generated  be  uniform   on    the    interval    (0,1),    and    that    the   sequence 
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be  non-repeating  for  a  sufficiently  long  period.   Other,  more 
involved  tests  also  exist.   This  method  has  the  advantage  that  it 
is  reproducible   (so  that  the  same  probabilities  may  be  regene- 
rated at  will  as  a  check  or  to  show  the  effect  of  different  initial 
conditions  on  the  final  outcome)  and  fast  (this  system  can  work 
"in  line") .   The  term  "congr uential"  refers  to  a  congruence  re- 
lation in  number  theory  and  modulo  arithmetic  (which  deals  only 
with  remainders)  Ll4]. 

Essentially,  the  method  relied  on  the  fact  that  the 
computer  (specifically,  the  IBM  36O/67)  had  a  fixed  storage  length 
of  four  bytes  (32  binary  digits  or  bits)  for  integer  constants. 
This  meant  that  the  highest  integer  the  computer  could  store  was 

2^^  -  1  =  2  lZf7  483  647 

a  ten  digit  number  (allowing  the  leading  bit  for  sign  plus  the 

31 
integer  "0",  there  were  2   -1  locations/combinations  for  other 

integers  [3].   If  one  multiplied  together  two  integers,  each  with 

five  or  six  digits,  there  was  generally  an  overflow  of  the  machine's 

storage  capacity.   FORTRAN  did  not  recognize  this  as  an  error,  and 

gave  only  the  last  eight,  nine  or  ten  digits;  that  is,  the  leading 

(and  generally  most  significant)  digits  were  truncated,  and  the 

remainder  was  stored  as  an  integer.   Essentially,  it  was  this 

integer  remainder  that  was  the  random  number.   To  make  it  useful,  one 

must  convert  it  to  a  floating  point  number  on  the  interval (0,1),  This 


was  done  by  making  the  integer  a  floating  point  number  and  dividing 

32 
by  2   ,  the  modulus  of  the  IBM  36O/67.   Since  about  half  the  integei 

remainders  (on  the  average)  were  negative  (due  to  a  negative  sign 


52 


in  the  sign  bit  of  the  integer  remainder),  one  actually  had  numbers 
on  the  interval  (-0.5,  +  0.5),  so  one  must  add  0.5  to  obtain  num- 
bers on  the  interval  (0,1)  fll]. 

In  practice,  a  total  of  four  FORTRAN  statements  were 
required  to  obtain  the  first  random  number.   They  were,  for 
example: 

9  KRAN    =  16807 
IRAN    =  87345 

10  IRAN    =  IRAN^KRAN 

11  RANDOM  =  0.5  +  FLOAT(IRAN)*2,328306  E-10. 

-32 
The  number  in  exponential  form  in  statement  11  is  2     written  as 

a  decimal  number .   Whenever  one  needed  a  new  random  number ,  one 
simply  used  statements  10  and  11.   Note  that  each  random  number 
"RANDOM"  was  generated  from  the  remainder  of  the  previous  integer 
remainder  "IRAN".   This  system  was  fast  since  it  worked  "in  line", 
and  computer  literature  indicates  that  this  generator  can  produce 
about  7,200  random  numbers  a  second  Cll]  .   Since  this  was  a  pseudo- 
random generator,  its  period  was  finite.   However,  the  theoretical 

29  r   T 

maximum  period  was  about  2    or  about  537  million  numbers  Ll4J. 

This  type  pseudo-random  generator  is  regularly  used  to  generate  a 

million  or  more  random  numbers.   Statistical  tests  performed  on 

the  numbers  used  and  observations  of  the  random  numbers  produced 

indicated  that  the  pseudo-random  numbers  generated  were  good  enough 

for  the  practical  testing  purposes  for  which  they  were  used  here. 

5.   Computer  Output 

The  problem  of  how  to  present  microstate  information  in 

easily  interpreted  form  was  solved  by  having  the  computer  "draw" 

pictures  of  microcrystal  planes,  called  "microstate  pictures"; 
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that  is,  the  computer  printed  arrays  of  zeroes  and  ones  to  indi- 
cate vacancies  or  occupied  sites,  and  each  digit  was  in  its  correct 
X  and  y  coordinate  position  for  the  site  it  represented.   See 
Figure  23.   The  program  contained  the  capability  of  printing  out 
microstate  pictures  for  IX  =  10 ,  IX  =  20  and  IX  =  26  (with  two 
planes  side  by  side),  and  for  IX  =  any  even  integer  through  60. 
The  only  restriction  on  lY  and  IZ  was  that  they  must  be  even 
integers.   The  program  also  contained  the  capability  of  choosing 
the  appropriate  output  segment  to  use. 

Variable  numbers  of  microstates  generated  were  allowed, 
and  this  was  provided  via  input  parameter  "NUMRUN"  -  the  number  of 
microstates  generated.   The  remainder  of  the  output  segment  was 
"nice  to  have"  information  that  was  printed  out  to  avoid  looking 
through  the  program. 

For  cases  in  which  large  numbers  of  microstates  were 
generated,  it  was  desired  to  have  only  summary  information  or 
microstate  pictures  every  mth  microstate,  and  have  all  microstate 
pictures  printed  after  a  certain  number  of  microstates  were  gene- 
rated.  This  capability  was  provided  by  input  parameters  "MSUM", 
"MPIX"  and  "MCRIT".   MSUM  was  the  index  indicating  that  a  short 
summary  of  the  current  microstate  was  printed  after  each  MSUM 
microstates  were  generated;  for  example,  if  MSUM  =  5,  a  summary 
was  printed  every  fifth  microstate.   The  summary  was  merely  the 
number  of  sites  occupied  in  each  z-layer.   MPIX  was  the  index 
indicating  that  a  microstate  picture  was  printed  after  each  MPIX 
microstates  were  generated;  for  example,  if  MPIX  =  10,  a  microstate 
picture  was  printed  of  every  tenth  microstate.   MCRIT  was  the  index 
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indicating   how  many   micr estates  were   generated  prior    to  printing    out 
pictures    of   every   microstate;    for    example,    if  MCRIT    =   250,    every 
microstate    from   250    to  NUMRUN   was   printed  as   a   picture. 
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APPENDIX  B 

LATTICE  LOADING  STATEMENTS 
The  different  initial  microstates  were  generated  by  sets  of 
one  or  more  program  statements  inserted  just  prior  to  statement 
45  in  the  lattice  generator.   The  statements  used  for  the  various 
lattices  generated  follow. 

1.  Perfect    Surface 

IF   (M-LL/2)      45,45,46 

This    statement    fills    the  bottom   half  of   the  active   lattice 
volume   with   atoms    and  leaves    the   top   half   empty. 

2.  Perfect    Surface   Plus    Randomly   Placed  Extra    Atoms 

IF   (M-LL/2)      45,45,49 
49  IF  (M-LL/2  -  ix*iY/2)   48,48,46 

48  IRAN  =  IRAN  *  KRAN 

RAND  =  FLOAT  (IRAN)  *    2.328306  E-10 
IF  (RAND)   45,45,46 

This  set  of  statements  fills  the  bottom  half  of  the  active 

lattice  with  atoms,  places  atoms  randomly  on  the  next  plane,  and 

leaves  all  other  active  lattice  planes  empty.   Sample  microstates 

achieved  with  this  set  of  statements  are  sketched  as  Figures  3a 

and  4a . 

3.  Perfect    Surface   Plus    Half  Plane   with   Edge   Parallel    to 
X-Axis 

IF    (M    -    (LL/2    +    IX  *   IY/4))    45,45,46 

This  statement  fills  the  bottom  half  of  the  active  lattice 
with  atoms,  places  a  half  plane  of  atoms  as  a  monatomic  step  in 
the  next  plane,  and  leaves  all  other  active  lattice  planes  empty. 
The  initial  microstate  achieved  is  sketched  as  Figure  5a. 
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4.  Pyramids 

lA  =  2 

IB  =  2 

IF  (LZ)  20,45,20 
20  IXP  =  LZ*(IA+1) 

IXM  =  IX-l-IXP 

IF  (LX-IXP)  46,25,25 
25  IF  (LX-IXM)  30,30,46 
30  lYP  =  LZ*(IB+1) 

lYM  =  lY-l-IYP 

IF  (LY-IYP)  46,35,35 
35  IF  (LY-iYM)  45,45,46 

This  set  of  statements  loads  a  pyramid  with  steps  that  move 

in  three  atomic  layers  at  a  time  -  see  Figure  6a.   To  achieve  a 

pyramid  with  steps  that  move  in  two  (one)  atomic  layers  (layer)  at 

a  time,  change  both  lA  and  IB  to  one  (zero)  -  see  Figure  7a. 

5.  Pyramid  holes 

lA  =  2 

IB  =  2 

IF  (LZ)  20,45,20 
20  IXU  =  (IX/2)  +  LZ*(IA+1)-1 

IXL  =  (IX/2)  -  LZ*(IA+1) 

IF  (LX-IXL)  45,25,25 
25  IF  (LX-IXU)  30,30,45 
30  lYU  =  (IY/2)  +  LZ*(IB+1)-1 

lYL  =  (IY/2)  -  LZ*(IB+1) 

IF  (LY-IYL)  45,35,35 

35  IF  (LY-IYU)  46,46,45 

This  set  of  statements  loads  a  pyramid  hole  with  steps  that 
are  three  atomic  layers  wide  -  see  Figure  8a.   To  achieve  a  pyramid 
hole  with  steps  two  (one)  atomic  layers  (layer)  wide,  change  both 
lA  and  IB  to  one  (zero)  -  see  Figure  9a. 

6.  Ridges 

To  achieve  a  ridge  parallel  to  the  y-axis,  use  the  same 
steps  as  for  pyramids,  but  use  the  statement 
IB  =  -1 
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step   width    is    determined  by   parameter    lA;    for : 

lA   =   2,    steps   are    three  atomic    layers    wide    -   see  Figure    Ha; 

lA    =    1,    steps    are   two  atomic    layers    wide; 

lA    =   0,    steps    are    one   atomic    layer    wide. 

To  achieve  a   ridge  parallel    to   the  x-axis,    again    use    the 
same   steps   as    for   pyramids,    but    use    the   statement 

lA   =    -1. 
Step   width    is    now   determined  by  parameter    IB;    for 

IB   =   2,    steps    are   three   atomic    layers    wide    -   see    Figure    10a; 

IB   =    1,    steps    are   two  atomic    layers    wide; 

IB   =   0,    steps    are    one   atomic   layer    wide. 

7.  valleys   Parallel    to   the   X-Axis 

IB   =    2 

IF  (LZ)   30,45,30 

30    lYU    =    (IY/2)    +   LZ*(IB+1)-1 
lYL    =    (IY/2)     -    LZ*(IB+1) 
IF    (LY-IYL)    45,35,35 

35  IF  (LY-IYU)   46,46,45 

This    set    of   statements    loads   a    valley   parallel    to    the   x-axis 
with   a   step   width   of    three   atomic    layers    -   see    Figure    12    a .      To 
achieve    such   a    valley    with   a    step   width    of   two    (one)    atomic    layers 
(layer),    change    IB   to    one    (zero)    -    see   Figure    13a. 

8.  Valleys    Parallel    to    the   Y-Axis 

lA   =   2 

IF    (LZ)    20,45,20 
20    IXU    =    (IX/2)    +    LZ*(IA+1)-1 
IXL    =    (IX/2)     -    LZ*(IA+1) 
IF    (LX-IXL)    45,25,25 

25  IF  (LX-ixu)   46,46,45 

This    set    of   statements    loads    a    valley   parallel    to   the  y-axis 
with   a   step   width   of   three  atomic    layers    -   see    Figure    l4a .      To 
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achieve  such  a  valley  with  a  step  width  of  two  (one)  atomic  layers 
(layer),  change  lA  to  one  (zero)  -  see  Figure  15a. 

9.  (1  -1  1)  Plane 

IF  (LX+LZ-LY)  45,45,46 

This    statement    loads    a    (1    -1    1)    plane    -   see   Figure    l6a. 

The  plane    is    truncated    if   z  is    smaller    than   either    x  or 

^  max  max 

y 

'max 

10.  (-1    1    1)    Plane 

IF   (LY+LZ-LX)    45,45,46 

This  statement  loads  a  (-1  1  1)  plane  -  see  Figure  17a. 

This  plane  is  truncated  if  z     is  smaller  than  either  x     or 

max  max 
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APPENDIX   C 

GLOSSARY  OF  COMPUTER  SYMBOLS 

Note:   This  Appendix  also  contains  the  definitions  of  several  terms 

peculiar  to  this  simulation. 

A(J,I):   A  4  X  5  array  coupled  between  the  MAIN  program  and  Sub- 
routine CROSYM;  it  carries  information  for  the  computation 
of  the  coefficients  of  the  cubic  potential. 

ALPHA:    Constant  read  in  as  data  used  in  computing  parameters  for 
the  truncated  Morse  potential. 

BOLTZ:    The  Boltzmann  constant,  equal  to  8.6171  x  lo"   eV/°K. 

CARRY  OUT  A  TRANSITION:   The  simulation  procedure  used  to  model  a 
"jump";  the  occupation  index  of  the  original  site  is  set 
equal  to  zero,  and  the  occupation  index  of  the  final  site 
(which  is  an  NNl  of  the  original  site)  is  set  equal  to  one. 

CFO:      A  constant  computed  by  the  program  for  use  in  the  force 
equation  resulting  from  the  cubic  potential. 

CFl :      A  constant  computed  by  the  program  for  use  in  the  force 
equation  resulting  from  the  cubic  potential. 

CF2 :      A  constant  computed  by  the  program  for  use  in  the  force 
equation  resulting  from  the  cubic  potential. 

CGBl :     A  constant  computed  by  the  program  for  use  in  the  truncated 
Morse  potential. 

CGB2:     A  constant  calculated  by  the  program  for  use  in  the  truncated 
Morse  potential. 
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CGDl:     A  constant  calculated  by  the  program  for  use  in  the  trun- 
cated Morse  potential. 

CGD2:     A  constant  calculated  by  the  program  for  use  in  the  trun- 
cated Morse  potential. 

CGFl :     A  constant  calculated  by  the  program  for  use  in  the  force 
equation  resulting  from  the  truncated  Morse  potential. 

03F2 :     A  constant  calculated  by  the  program  for  use  in  the  force 
equation  resulting  from  the  truncated  Morse  potential. 

CCX^A:     Labeled  COMMON  storage. 

CPO:      A  constant  computed  by  Subroutine  CROSYM  for  use  in  the 
cubic  potential. 

CPl:      A  constant  computed  by  Subroutine  CROSYM  for  use  in  the 
cubic  potential. 

CP2 :      A  constant  computed  by  Subroutine  CROSYM  for  use  in  the 
cubic  potential. 

CP3:      A  constant  computed  by  Subroutine  CROSYM  for  use  in  the 
cubic  potential. 

CROSYM:   A  Subroutine  developed  at  Lockheed  to  solve  several 

simultaneous  linear  equations  by  the  "Method  of  Crout". 

CVD:      "Distance  conversion  factor",  CVR  x  10~  to  convert  lattice 
units  to  meters . 

CVDE:     CVD/CVE,  a  ratio  to  avoid  repeated  division;  the  reciprocal 

of  CVED. 

-19 
CVE:  "Energy   conversion   factor",    1.6   x    10  to   convert     elec- 

tronvolts    to   joules. 

CVED:  CVE/CVD,    a    ratio   used   to  avoid  Repeated  division;    the 

reciprocal    of   CVEE. 
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-27 
CVM:      ''Mass  conversion  factor",  I.672  x  10     to  convert 

atomics  mass  units  to  kilograms 

CVR :      "Radial  distance  conversion  factor",  the  lattice  unit  in 

angstroms,  used  to  convert  lattice  units  to  angstrom  units. 

DCON:     Constant  read/in  as  data,  and  used  in  computing  parameters 
for  the  truncated  Morse  potential. 

DE(J,I):  "Difference  in  energy",  an  array  containing  the  energy 
difference  between  an  atom  at  site  I  and  its  Jth  NNl 
position,  whether  occupied  or  vacant  (if  occupied,  the 
particular  array  value  is  made  arbitrarily  large  to  avoid 
matter  destruction);  the  values  on  this  array  are  used 
in  determining  transition  probabilities;  this  array  must 
be  dimensioned  at  least  12  x  LL. 

DIS:      SQRT(DIST2(I) ) 

DIST(I):  "Distance";  a  vector  array  containing  the  value  of  the 
distance  (in  lattice  units)  between  Ith  NN's. 

DIST2(I):  "Distance  squared";  a  vector  array  containing  the  value 

of  the  distance  squared  (in  lattice  units  squared)  between 
Ith  NN's. 

EMIN:     "Minimum  energy";  the  value  at  the  lowest  DE(J,I),  found 

by  searching,  when  a  jump  is  to  be  made  to  a  site  at  lower 
potential;  o£  the  current  value  of  DE(J,I),  found  by 
searching  and  working  from  lowest  to  highest  energy  dif- 
ferences, used  when  attempting  a  jump  to  a  site  of  higher 
potential. 
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EMINXT:   "Next  minimum  energy";  the  value  of  the  next  lower  DE(J,I), 
found  by  searching,  used  when  attempting  to  jump  to  a  site 
of  higher  potential  after  all  attempts  at  jumping  have 
failed  for  lower  positive  energy  differences. 

EPSLON:   "Epsilon";  a  constant  used  to  determine  if  two  DE(J,I)'s 
were  effectively  equal;  its  value  is  1.0  x  10    eV. 

EQTEMP:   "Equivalent  temperature";  the  effective  lattice  temperature 
achieved  as  a  result  of  varying  PFIV;  EQTEMP  =0.5 
*  TEMP/PFIV,  and  this  relation  is  plotted  in  Figure  22. 

EV:       "Electron  volt". 

EXA:      A  constant  read  in  as  data  and  used  in  calculation  of  the 
Gibson  II  potential. 

EXB:      A  constant  read  in  as  data  and  used  in  calculation  of  the 
Gibson  II  potential. 

FLAT:     Constant  read  in  as  data  and  equal  to  2.0*CVR. 

FRCAND(X):  "Anderraan  force";  a  function  defined  in  the  MAIN  program 
for  calculation  of  the  force  arising  due  to  the  truncated 
Morse  potential. 

FRCGIB(X) :  "Gibson  force";  a  function  defined  in  the  MAIN  program 
for  calculation  of  the  force  arising  due  to  the  Gibson  II 
potential . 

FTEST(I,J):  A  12  x  12  array  used  in  choosing  the  site  to  which  an 
atom  jumps  if  more  than  one  of  the  NNl  sites  of  an  atom 
under  consideration  have  the  same  energy;  for  different 
crystal  orientations  the  number  of  NNl's  may  change,  and 
therefore  the  dimensions  of  array  FTEST  must  also  change. 
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FXA: 


I: 

lA: 


IB: 


I COUNT: 


IFUL; 


IH2: 


IHEAD: 


IHNOl : 


IHTARG : 


IHTPOr ; 


A  constant  computed  by  the  program  for  use  in  the  force 
equation  resulting  from  the  Gibson  II  potential. 
A  DO  loop  index. 

A  constant  used  by  the  lattice  loader  to  determine  the 
step  width  of  the  initial  micr estate  for  several  different 
configurations;  see  Appendix  B. 

A  constant  used  by  the  lattice  loader  to  determine  the 
step  width  of  the  initial  microstate  for  several  different 
configurations;  see  Appendix  B. 

A  counting  index  which  is  the  number  of  the  microstate 
being  generated  by  the  program.  i 

"Full";  the  number  of  a  lattice  site  in  the  perfect  cry- 
stalline substrate  on  which  the  active  lattice  volume  is 
placed;  all  lattice  sites  in  this  substrate  are  always 
full  and  each  is  represented  by  site  number  IFUL. 
Alphanumeric  array  read  in  as  data  for  truncated  Morse 
potential  parameters. 

Alphanumeric  array  read  in  as  data  to  provide  an  output 
heading  for  the  particular  simulation  attempted. 
Alphanumeric  array  read  in  as  data  to  provide  a  program 
flag  warning  that  during  a  calculation  an  atom  was  found 
that  had  no  NNl's  (this  is  unphysical  in  this  simulation); 
the  message  merely  states  "NO  NNl's". 

Alphanumeric  array  read  in  as  data  and  used  to  define  the 
particular  metal  simulated. 

Alphanumeric  array  read  in  as  data  to  provide  a  title  for 
Gibson  II  potential  parameters. 
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II:  A  DO  loop  index. 
IJ:  A  DO  loop  index. 
IK:  A  DO  loop  index. 
IRAN:     The  present  (or  initial)  value  of  the  constant  used  in 

computing  pseudo-random  numbers. 
IRANI:    ^'Initial  IRAN";  defined  to  retain  a  knowledge  of  the 

original  random  number  seed  used  in  a  simulation. 
ISX:      IX/2,  a  shift  index  indicating  a  shift  of  one  row  of  atoms. 
ISXM:     ISX-1,  a  shift  index  indicating  a  shift  of  one  atom  less 

than  a  full  row. 
ISXP:     ISX+1,  a  shift  index  indicating  a  shift  of  one  atom  more 

than  a  full  row. 
ISZ:      ISX*IY,  a  shift  index  indicating  a  shift  of  one  plane  of 

atoms . 
ISZM:     ISZ-1,  a  shift  index  indicating  a  shift  of  one  atom  less 

than  a  full  plane. 
ISZP:     ISZ+1 ,  a  shift  index  indicating  a  shift  of  one  atom  more 

than  a  full  plane. 
IT:       An  index  used  in  computing  ITT. 
ITEST:    A  DO  loop  index;  also,  a  constant  used  to  determine  which 

site  may  receive  an  atom  when  carrying  out  a  transition 

in  the  case  of  several  available  sites  with  the  same 

energy. 
ITT:      IT+JT+KT,  the  index  which  determines  when  or  not  a  particular 

coordinate  position  should  be  labeled  as  a  valid  lattice 

site  by  the  lattice  generator. 
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IVAC:  "Vacant";    the   number    of   a    lattice    site    in    the   void  above 

the  active    lattice   volume;    all    sites    in   this    void  are 

always    vacant   and   each    is    represented  by    site   number    IVAC. 
IX:  An    even     integer  read   in   as    data    indicating    the    number    of 

active    lattice    volume   planes    generated   in    the   x-direction. 
IXL:  A   variable   used  as    a    shifter    when   generating  pyramid   holes 

and   valleys;    see   Appendix   B. 
IXM:  A   variable   used  as    a    shifter    when    generating   pyramids    and 

ridges;    see   Appendix    B. 
IXP:  A   variable   used  as   a    shifter    when    generating   pyramids   and 

ridges;    see   Appendix    B. 
IXU:  A   variable  used   as    a    shifter    when   generating  pyramid  holes 

and   valleys;    see   Appendix    B. 
lY:  An    even    integer    read    in   as    data    indicating    the   number    of 

active    volume   planes    generated    in    the   y-direction. 
lYL:  A   variable   used  as   a    shifter    when    generating . pyramid  holes 

and  valleys;    see  Appendix  B. 
lYM:  A   variable   used  as    a    shifter    when   generating   pyramids   and 

ridges;    see   Appendix   B. 
lYP:  A   variable   used  as    a    shifter    when    generating  pyramids   and 

ridges;    see  Appendix    B. 
lYU:  A   variable   used  as    a    shifter    when    generating   pyramid  holes 

and   valleys;    see   Appendix    B. 
IZ:  An    even    integer    read   in   as    data    indicating    the   number    of 

active   lattice   volume  planes    generated   in    the  z-direction. 
J:  A   DO   loop    index. 
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jj:       lY-J,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
JJM:      JJ-1,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
JLOE(I):  "Low  energy  index";  a  vector  array  used  to  record  the  site 

number  of  the  lowest  energy  ^fNl  site  of  an  atom  attempting 

to  make  a  transition;  this  vector  must  have  dimension  of 

at  least  12. 
JMAX:     J-lj  an  index  indicating  the  number  of  NNl  sites  of  lower 

energy  for  an  atom  attempting  to  make  a  transition. 
JPl:      J+1,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
JT:       An  index  used  in  computing  ITT. 
JUMP:     The  physical  process  of  an  atom  leaving  its  present 

location  and  moving  into  a  vacant  NNl  site. 
KA:       K*ISZ,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
KB:       KA-J*ISZ+1,  a  constant  used  in  the  output  segment  to  draw 

and  label  microstate  pictures. 
KC:       KA-JP1*ISX+1,  a  constant  used  in  the  output  segment  to  draw 

and  label  microstate  pictures. 
KD:       KB+ISXM,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
KB:       KC+ISXM,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
KF:       KP1*ISZ,  a  constant  used  in  the  output  segment  to  draw  and 

label  microstate  pictures. 
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KG:      KF-J*ISX+1,  a  constant  used  in  the  output  segment  to  draw 

and  label  microstate  pictures. 
KH:       KF-JP1*ISX+1 ,  a  constant  used  in  the  output  segment  to 

draw  and  label  microstate  pictures. 
KI:      KG+ISXM,  a  constant  used  in  the  output  segment  to  draw 

and  label  microstate  pictures. 
KJ:       KH+ISXM,  a  constant  used  in  the  output  segment  to  draw 

and  label  microstate  pictures. 
KKOUNT:   An  index  used  to  avoid  recalculating  potential  coefficients 

when  more  than  one  crystal  simulation  is  attempted. 
KMl :      K-l,an  index  used  in  printing  out  microstate  summary  infor- 
mation; also,  a  constant  used  to  draw  and  label  microstate 

pictures. 
KNEW:     "New  K";  site  to  which  a  transition  has  just  been  carried 

out. 
KOUNTM:   K0UNT2-1,  an  index  used  in  deciding  which  site  to  jump 

to  when  there  is  more  than  one  NNl  site  of  lowest  energy 

to  choose  from. 
KOUNTl :   A  counting  index  used  to  determine  the  number  of  NNl  sites 

that  have  tower  or  higher  energy  with  respect  to  the  site 

under  consideration. 
K0UNT2 :   A  counting  index  used  to  determine  the  multiplicity  of 

NNl  sites  with  lowest  energy. 
K0UNT3:   A  counting  index  used  to  determine  the  multiplicity  of 

NNl  sites  with  higher  energy. 
K0UNT4:   A  counting  index  used  to  determine  whether  or  not  all  higher 

energy  sites  have  been  tested  for  the  possibility  of  making 

a  transition. 
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K0UNT5: 
K0UNT6 : 
KPl: 

KP2: 
KP3: 

KRAN: 

KT: 

LATTICE 

LDX: 

LDY: 

LDZ: 

LL: 

L.U  : 
LX: 


An    index  controlling  program  flow  when  attempting   to 
determine    if  all    sites    of  higher    energy   have  been    tested. 
A   counting    index   indicating    the   multiplicity    of   NNl    sites 
with   next    higher    energy. 

K+1,    an    index   used   in   printing   out   microstate   summary 
information;    also,    a    constant    used   to   draw  and   label 
microstate  pictures. 

K+2,    an    index   used   in  printing   out    microstate   summary 
information . 

K+3,    an    index  used    in   printing   out   microstate   summary 
information. 

A   constant    used   in    computing  pseudo-random   numbers. 
An    index   used   in    computing    ITT. 
UNIT.(L.U.):    The    distance   between    two  NNl's    along   a    single 
coordinate    direction. 

A   shift    index  used   to  assign    lattice   site    coordinate 
positions    in    the   x-direction. 

A   shift    index   used   to   assign    lattice   site   coordinate 
positions    in    the   y-direction. 

A   shift    index   used   to  assign    lattice    site   coordinate 
positions    in    the   z-direction. 

"Lattice    length";    the   number    of  available   sites    in    the 
active    lattice    volume,    equal    to    IX*IY*IZ/2, 
See   LATTICE   UNIT. 

The   x-coordinate  position   assigned   to  a   particular    site 
during    lattice   generation. 
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LY:      The  y-coordinate  position  assigned  to  a  particular  site 
during  lattice  generation. 

LZ:       The  z-coordinate  position  assigned  to  a  particular  site 
during  lattice  generation. 

M:        A  counting  index  used  to  assign  a  label  to  lattice  sites; 
also,  an  index  (=  MCRO)  used  in  Subroutine  CROSYM . 

MCRIT:    A  constant  read  in  as  data;  all  microstates  generated  from 
microstate  number  MCRIT  to  NUMRUN  are  printed  out  as 
microstate  pictures. 

MCRO:     A  constant  read  in  as  data,  indicating  the  number  of 
unknowns  to  be  solved  for  by  Subroutine  CROSYM. 

MICROSTATE:  The  surface  configuration  either  initialized  in  the 
program  or  generated  by  it. 

MICROSTATE  PICTURE:  An  array  of  occupation  indices  printed  out  in 
their  correct  coordinate  positions  to  give  a  physical 
picture  of  the  surface  configuration  (see  Figure  23). 

MICROSTATE  SUMMARY:  A  listing  of  the  microstate  number  and  the 
number  of  atoms  in  each  z -plane. 

MPIX:     A  constant  read  in  as  data  indicating  that  between  micro- 
state  zero  and  microstate  MCRIT,  a  microstate  picture  is 
printed  out  every  MPIXth  microstate. 

MSUM:     A  constant  read  in  as  data,  indicating  that  between  micro- 
state  zero  and  microstate  MCRIT,  a  microstate  summary  is 
printed  every  MSUMth  microstate  (except  when  a  microstate 
picture  is  printed). 
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MULT(I):  ^'Multiplicity";  a  vector  array  indicating  the  numbers  of 
NNl  sites  which  have  the  same  energy  difference  with 
respect  to  the  atom  of  interest;  this  vector  must  have 
'   dimension  at  least  as  great  as  the  number  of  NNl's  for  a 
lattice  location  (in  this  simulation,  the  proper  number 
was  12)  . 

NA:       An  index  indicating  which  site  an  atom  jumped  to  when 

several  NNl  sites  of  the  same  energy  were  located;  the 
value  of  NA  was  chosen  stochastically. 

NBRFOR(J,K):  "Neighbor  four";  an  array  containing  the  site  numbers 
of  the  12  NN4's  of  the  Ith  lattice  site;  since  this  simu- 
lation did  not  use  NN4*s  the  array  was  dimensioned  1x1, 
but  if  Nn4's  are  included,  it  must  have  dimensions  of  at 
least  12  X  LL. 

NBRONE(J,I):  "Neighbor  one";  an  array  containing  the  site  numbers 
of  the  12  NNl's  of  the  Ith  lattice  site;  this  array  must 
have  dimensions  of  at  least  12  x  LL. 

NBRTHR(J,I):  "Neighbor  three";  an  array  containing  the  site  numbers 
of  the  24  NN3's  of  the  Ith  lattice  site;  since  this  simu- 
lation did  not  use  NN3's,  the  array  was  dimensioned  1x1, 
but  if  NN3's  are  included,  it  must  have  dimensions  of  at 
least  24  X  LL. 

NBRTWO(J,I):  "Neighbor  two";  an  array  containing  the  site  numbers 
of  the  six  NN2's  of  the  Ith  lattice  site;  this  array  must 
have  dimensions  of  at  least  6  x  LL. 
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NCRYST:   An  index  read  in  as  data  to  indicate  the  number  of  crystal 
simulations  to  be  run  using  a  single  input  deck;  if 
NCRYST  =  0,  the  simulation  with  which  this  value  is  read 
in  is  the  last  (or  the  only)  simulation  run;  if  NCRYST  7^  0, 
there  is  at  least  one  simulation  following  the  one  con- 
taining a  non-zero  value  of  NCRYST. 

NN:       "Nearest  neighbor". 

NNl:      "First  nearest  neighbor". 

NN2 :      "Second  nearest  neighbor". 

NNDIS2(I):  "Nearest  neighbor  distance  squared";  the  square  of  the 
distance  to  the  Ith  NN;  this  vector  must  have  dimension 
at  least  as  great  at  NUMNN. 

NNFX30C(I):  "Number  of  Nn4  sites  occupied",  defined  with  respect  to 
site  I;  this  vector  must  have  dimension  at  least  as  great 
as  LL. 

NNONOC(I):  "Number  of  NNl  sites  occupied",  defined  with  respect 
to  site  I;  this  vector  must  have  dimension  at  least  as 
great  as  LL. 

NNTHOC(I):  "Number  of  NN3  sites  occupied",  defined  with  respect  to 
site  I;  this  vector  must  have  dimension  at  least  as  great 
as  LL. 

NNTWOC(I):  "Number  of  NN2  sites  occupied",  defined  with  respect  to 
site  I;  this  vector  must  have  dimension  at  least  as  great 
as  LL. 

NOCC(I):  Vector  array  containing  the  occupation  index  of  each 

lattice  site;  this  vector  must  have  dimension  at  least 
as  great  as  LL+3. 
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NPLANE:    "Plane   number";    NZ(I)+1,    a    constant    used   to   define  plane 

numbers    greater    than    zero   for    use   as    subscripts. 
NUMATM :     "Number    of   atoms";    the    number    of   atoms    actually  placed    in 

the   active    lattice    volume   by    the    lattice    generator. 
NUMNN:     "Number    of  nearest    neighbors",    considered    in   a   particular 

simulation. 
NUMRUN:       "Number    of  runs";    a    constant    read   in   as    data    indicating 

the   number    of   microstates    to   be   generated  by    the    simu- 
lation. 
NUMSIT(I):     "Number    of   sites"   occupied    in    the    Ith  plane;    here,    I    is 

defined  by    NPLANE    in    order    to   have   non-zero   subscripts. 
NX(I):         Vector    array    containing    the   values    of   the   x-coordinate 

for  each  lattice  site;  this  vector  must  have  dimension 

at  least  as  great  as  LL. 
NY(I):     Vector  array  containing  the  values  of  the  y-coordinate 

for  each  lattice  site;  this  vector  must  have  dimension 

at  least  as  great  as  LL. 
NZ(I):     Vector  array  containing  the  values  of  the  z-coordinate 

for  each  lattice  site;  this  vector  must  have  dimension 

at  least  as  great  as  LL. 
OCCUPATION  INDEX:  A  number  indicating  whether  or  not  a  site  is 

occupied;  an  occupation  index  of  zero  implies  that  the 

site  is  vacant,  and  an  occupation  index  of  one  implies 

that  the  site  is  occupied;  no  other  numbers  are  allowed 

as  an  occupation  index. 
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PE(I):    "Potential  energy";  a  vector  array  containing  the  value  of 

the  potential  energy  at  each  lattice  site  due  to  nearby 

occupied  sites;  this  vector  must  have  dimension  at  least 

as  great  as  LL . 
PENN(I):  "Potential  energy  due  to  nearest  neighbor";  the  potential 

energy  contribution  due  to  the  Ith  NN. 
PFIV:     "Point  five";  an  energy  distribution  factor  employed  in 

calculating  site  potential  energy  and  equivalent  lattice 

temperature;  the  correct  average  value  of  PFIV  is  0.5. 
POTF(X):  A  function  defined  in  the  MAIN  program  for  calculation 

of  the  Gibson  II  potential. 
P0T2F(X):  A  function  defined  in  the  MAIN  program  for  calculation 

of  the  cubic  potential. 
P0T3F(X):  A  function  defined  in  the  MAIN  program  for  calculation  of 

the  truncated  Morse  potential. 
RAND:     FLOAT( IRAN)*2 .328306E-IO,  a  pseudo-random  number  on  the 

interval  (-0.5,  0.5)  used  in  placing  atoms  randomly  on  a 

perfect  surface;  see  Appendix  B. 
RANDOM:   0 . 5+FL0AT( IRAN) *2 . 328306E-IO ,  a  pseudo-random  number  on 

the  interval  (0,1)  used  for  comparison  with  probabilities 

when  a  decision  based  on  a  probability  was  made. 
RE:       A  constant  read  in  as  data  used  in  computing  truncated 

Morse  potential  parameters. 
ROEA:     A  constant  read  in  as  data  equal  to  the  distance  in 

lattice  units  at  which  the  Gibson  II  and  cubic  potentials 

match. 
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ROEB:     A  constant  read  in  as  data  equal  to  the  distance  in 

lattice  units  at  which  the  cubic  and  truncated  Morse 

potentials  match. 
ROEC:     A  constant  read  in  as  data  equal  to  the  distance  in 

lattice  units  at  which  the  truncated  Morse  potential  is 

truncated. 
R0EC2:    "ROEC  squared";  ROEC*ROEC. 
SUM:      A  variable  used  by  Subroutine  CROSYM . 
TAR:      Alphanumeric  array  read  in  as  data  and  not  used  in  this 

simulation. 
TAU:      BOLTZ*TEMP,  i.e,  the  "kT"  factor  used  in  computing  a 

Boltzmann  factor. 
TEMP:     ^'Temperature";  the  nominal  lattice  temperature;  see 

EQTEMP . 
TMAS:     "Mass";  a  constant  read  in  as  data  equal  to  the  mass  of 

a  lattice  atom  in  atomic  mass  units. 
TPOT:     "Total  potential";  the  sum  of  all  occupied  site  potentials 

due  to  interactions  with  nearest  neighbors  through  NUMNN. 
TPROB:    "Transition  probability",  defined  by  equations  A-7. 
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COMPUTER    PROGRAM 
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TABLE  I 


POTENTIAL  FUNCTIONS  AND  COEFFICIENTS 


^Gibson  li^)  =  exp(EXA-HEXB-R) 


V    .  (R)  =  R*(R^(R*CP3+CP2)+CP1)+CP0 

V^  „  (R)    =    exp(CGDl+CGBl*R)     -    exp(OGD2+CGB2*R ) 

trun  Morse^     /  f\  /  f\  i 


Type   Potential 

Coefficients 

Cutoff   Distance 
( lattice   units ) 

Gibson    II 

EXA      =    10.02407 
EXB      =    -9.19667 

'a  =  °-«3 

cubic 

cpo    =  587.6182 

'b  =  1-1° 

CPI     =  -1593.863 

CP2      =    1450.286 

CP3     =   -442.2820 

truncated  Morse 

CGDl    =    6.649136 
CGD2    =    3.651313 
CGBl    =    -5.077808 
CGB2    =    -2.538904 

'c  "  ^•'*° 
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Figure  1.   The  Active  Lattice  Volume. 


GENERATE  LATTICE 

load  initial  microstate 
provide  periodic  boundary  conditions 

' 

COMPUTE  COMPOSITE  POTENTIAL  FUNCTION 
calculate  potential  for  each  NN 
assign  each  site  a  potential  energy 

COMPUTE  TRANSITION  PROBABILITIES 
decide  whether  or  not  an  atom  jumps 
if  an  atom  does  jump,  patch  potential  of 
nearby  atoms 

NEW  MICROSTATE 
print  out  desired  information 
recycle  to  generate  another  microstate 

Figure  2.  Overall  Program  Flow 
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Figure   3.    Perfect    Surface   Plus    99   Randomly   Placed  Atoms. 
(a)    Initial   raicrostate    (lines    bordering   atoms 
and  atom  clusters    are  parallel   and  perpendicular 
to   the    (lOO)    direction). 


[b)    Microstate    l40    (lines   bordering   atom   clusters    are 
parallel   and  perpendicular    to   the    (lio)    direction] 
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Figure  4.  Perfect  Surface  Plus  104  Randomly  Placed  Atoms. 
(a)  Initial  microstate  (lines  bordering  atoms 
and  atom  clusters  are  parallel  and  perpendicular 
to  the  <100>  direction). 


(b)  Microstate  l40  (lines  bordering  atom  clusters  are 
parallel  and  perpendicular  to  the  ^lio)  direction). 
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Figure  5.  Perfect  Surface  Plus  Half  Plane  Monatomic  Step; 
(a)  initial  microstate 


Figure  5b.  Microstate  500  (lines  bordering  the  atom 
cluster  are  parallel  and  perpendicular 
to  the  <110>  direction) 
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Figure  6.  Pyramid  With  Steps  Three  Atomic  Layers  Wide; 
(a)  initial  microstate. 


Figure  6b.  Microstate  500  (lines  bordering  the  atom 
cluster  are  parallel  and  perpendicular  to 
the  (110>  direction). 
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Figure  7.  Pyramid  With  Steps  Two  Atomic  Layers  Wide. 
(a)  Initial  microstate. 


(b)  Microstate  500  (lines  bordering  atom  clusters  are 
parallel  and  perpendicular  to  the  (110>  direction);  the 
hole  in  the  z  =  1  plane  is  a  result  of  the  partial  fill-in 
and  reorientation  of  the  initial  grooves  in  the  z  =  1  plane; 
however,  the  vacancy  in  the  z  =  0  plane  (lower  left  corner) 
is  due  to  an  atom  jumping  out  of  a  perfect  surface  (the 
z  =  0  plane)  -  this  phenomenon  was  observed  infrequently. 
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Figure  8.  Pyramid  Hole  With  Steps  Three  Atomic  Layers 
Wide;  (a)  initial  microstate. 


Figure  8b.  Microstate  500  (lines  bordering  atoms  and 

atom  clusters  are  parallel  and  perpendicular 
to  the  <110>  direction). 
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Figure   9.    Pyramid  Hole   With  Steps    Two  Atomic  Layers   Wide, 
(a)    Initial   microstate. 


(b)    Microstate    500    (lines    bordering   clusters    of  atoms    are 
parallel   and  perpendicular    to   the   <110)    direction);    the 
holes    appearing    here  are   due    to   the  partial    fill-in   and 
reorientation    of    the   initial   holes. 
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/- 


Figure    10.    Ridge  Parallel    to   X-Axis   With    Steps   Three 
Atomic   Layers   Wide,    (a)    Initial   raicrostate, 


(b)   Microstate    500    (lines    bordering   atom    clusters    are 
parallel    and  perpendicular    to   the   (llO)    direction);    the 
hole   appearing    in    the   z    =    1   plane    is    due    to   the  partial 
fill-in   and   reorientation    of   the    initial    groove    in    the 
z    =    1   plane. 
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Figure   11.    Ridge  Parallel    to  Y-Axis    With  Steps   Three 

Atomic  Layers    Wide,    (a)    Initial   microstate. 


(b)    Microstate    500    (lines    bordering   atom   clusters    are 
parallel    and  perpendicular    to   the   (lio)    direction);    the 
hole   appearing    in    the   2=1   plane    is   a    result    of    the 
partial    fill-in   and   reorientation   of   the    initial    groove 
in   the    z    =    1  plane.       ii4 


(c)  Microstate  451,  using  a  different  random  number  seed 
then  used  to  achieve  the  microstate  depicted  in  Figure  lib 
(the  lines  bordering  atoms  and  atom  clusters  are  parallel 
and  perpendicular  to  the  (lio)  direction);  the  hole 
appearing  in  the  2=1  plane  is  a  result  of  the  partial 
fill-in  and  reorientation  of  the  initial  groove  in  the 
z  =  1  plane. 
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Figure    12.    Valley   Parallel    to  X-Axis   With   Steps   Three 
Atomic   Layers   Wide;    (a)    initial   microstate 


Figure    12b.    Microstate   500    (lines   bordering    the   atom 

clusters   are  parallel   and  perpendicular    to 
the    <110>    direction).      The   hole  appearing 
here    is    a   result    of    the  partial    fill-in   and 
reorientation    of   the    initial    groove    in    the 
z=l   plane. 
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Figure   13.    Valley   Parallel    to  X-Axis   With   Steps   Two 

Atomic  Layers   Wide;     (a)    initial   microstate. 


Figure   13b.    Microstate    500    (lines    bordering    the  atom 

clusters    are  parallel    and  perpendicular    to 
the   (110>    direction).      The   hole   appearing 
here    is    a    result    of    the  partial    fill-in 
and  reorientation    of   the    initial    groove    in 
the  2,-2   plane. 
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Figure    l4.    Valley   Parallel    to   Y-Axis    With   Steps    Three 
Atomic   Layers   Wide,     (a)    Initial   microstate, 


(b)    Microstate    500    (lines    bordering   atom   clusters    are 
parallel   and  perpendicular    to  the   <110>    direction);    the 
hole   appearing    in    the   z    =    1   plane    is    a    result      of   the 
partial    fill-in   and  reorientation   of    the    initial    groove 
in    the  z    =    1   plane. 
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Figure   15.    Valley   Parallel    to   Y-Axis    With    Steps   Two  Atomic 
Layers   Wide,    (a)    Initial   raicrostate. 


(b)    Microstate    500    (lines    bordering   atom   clusters    are 
parallel   and  perpendicular    to   the   (llO>    direction). 
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Figure  l6 .  (1  -1  1)  Plane,  Truncated,  (a 
micr estate , 


Initial 


(b)    Microstate    500. 
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Figure   17.    (-1    1    1)    Plane,    Truncated,    (a)    Initial 
micr estate. 


(b)    Microstate    500, 
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Figure   l8.    The   12   NNl   Positions 

(a)  In    the   z-plane   above 
the   atom; 

(b)  In    the   same   z-plane  as 
the   atom; 

(c)  In   the  z-plane   below 
the  atom. 
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Figure    19.    The  Six   NN2 
Positions . 

Sites    at   positions    1-4  are 
->"X         in   the   same   z-plane  as    the 
atom;    position    5    is    two 
z-planes    below    the   atom; 
position   6    is    two   z-planes 
above   the  atom. 
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Figure  20.  The  Potential  Functions 
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Figure  21.  The  Composite  Potential 
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Figure  22.  Equivalent  Lattice  Temperature  as  a 
Function  of  PFIV. 
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Figure  23.  Sample  Microstate  Picture.   Here,  IX  =  10 
and  lY  =  4. 
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Figure   24.    Microstate   310    of   an    Initially   Perfect    Surface 
at    10,000    °K.         (a)    Lines    bordering   atoms   and 
clusters    of  atoms    are  parallel   and  perpendicular 
to  the    (lOO)    direction. 


(b)    Lines    bordering   atoms    and   clusters    of   atoms   ar< 
parallel    and  perpendicular    to   the   (lio)    direction. 
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Figure   25.    Microstate    l80   of  an    Initially   Perfect    Surface 
at   2,000   °K    (lines    bordering   atom   clusters    are 
parallel   and  perpendicular    to   the   <llo) 
direction) .    Note   that    atoms    have    jumped    out    of 
the  perfect    surface  and   coalesced   on   the   surface, 


Figure 


Microstate    190   of  an    Initially   Perfect   Surface 
at    2,584    ^    (lines    bordering   atom  clusters    are 
parallel    and  perpendicular    to   the    <110) 
direction).    Note    that   a    hole   two  atomic    layers 
deep    has    been   burroughed  as    atoms    left    the 
perfect    surface  and  piled  up    near    it. 
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Figure  28.  Break  Repair  Sequence  For  Atom  Aggregates  in  a 
Single  Plane.  The  initial  microstate  was  Figure 
10a. 


Figure  28a.  Microstate  50,  z  =  2  plane  (75  atoms) 


Figure  28b.  Microstate  100,  z  =  2  plane  (70  atoms) 


Figure  28c.  Microstate  120,  z  =  2  plane  {66    atoms). 
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Figure  28d.  Microstate  130,  z  =  2  plane  {66   atoms). 


Figure  28e.  Microstate  l40,  z  =  2  plane  (64  atoms 


Figure  28f.  Microstate  150,  z  =  2  plane  (63  atoms).  The 
two  groups  of  atoms  are  held  together  only 
by  one  NN2  interaction. 
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Figure  28g.  Microstate  170,  2=2  plane  (6l  atoms). 


Figure  28h.  Microstate  l80,  z  =  2  plane  (6l  atoms) 


Figure  28io  Microstate  190,  2=2  plane  (60  atoms).  The 
two  groups  of  atoms  are  joined  only  by  one 
NN2  interaction. 
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Figure  28j.  Microstate  200,  2=2  plane  (59  atoms). 


Figure  28k.  Microstate  240,  z  =  2  plane  (57  atoms). 


Figure  28l.  Microstate  290,  7.-2   plane  (63  atoms  -  i.e. 
some  atoms  have  fallen  down  from  the  z  =  3 
plane) . 
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Figure  28ra.    Microstate   330,    z    =   2  plane    (6l   Atoms). 


Figure   28n.    Microstate    350,    z    =   2   plane    (57   atoms) 


Figure   28o.    Microstate   370,    2=2  plane    [55  atoms) 
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Figure   29.    Relative  Atomic  Positions    in  a    Lattice   Plane, 
Here,    IX   =  4  and   lY   =  6 . 
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13.    ABSTRACT 


A   simple  physical   model    of   a    copper    crystal   surface   was    developed.      Atoms 
were   considered   as    quasi-hard   spheres    which    occupied  perfect    lattice  positions 
A   computer    simulation,    based   on   energy    consideration    only,    using    the   Monte 
Carlo   method  was    developed,    tested  and   used   to   study    equilibrium   surface 
micr estates.      As    a    result    of   this    study,    four    conclusions    were    drawn: 

1.  This  model  holds  promise  for  further  investigation  of  real  crystal 
surface  phenomenon. 

2.  Minimum  energy  considerations  cause  atoms  to  align  themselves  pre- 
ferentially in  a  (no)  direction  on  the  surface  of  a  face-centered  cubic 
crystal . 

3.  Stepped   surface    configurations   are   fairly    stable,    but    isolated    "stub" 
atoms   and   vacancies    tend   to   coalesce   with    other    "stubs"   and   vacancies, 
respectively . 

4.  Random   motions    of   the    individual   atoms    cause   aggregates    of    atoms    to 
break  apart    and   recombine. 
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